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1.     INTRODUCTION 

This  report  deals  with  some  calibration  problems  of  monitoring  a 
multiple  array  underwater  tracking  range.  The  arrays  in  the  system  are  of  the 
short  baseline  type;  each  contains  four  sonar  transducers  placed  rigidly  at  the 
corners  of  a  cube  in  a  manner  that  describes  a  Cartesian  coordinate  system  in 
three  dimensions.  Figure  1  contains  a  diagram  showing  the  structure  for 
these  arrays  and  indicating  the  numerous  signals  they  may  receive.  It  is 
slightly  deceptive  in  that  real  ray  paths  are  not  straight  lines. 

An  array  receives  a  distinctive  signal  from  a  synchronously  timed  pinger 
attached  to  the  target  vehicle.  The  differentials  of  the  sound  wavefront's 
times  of  arrival  at  the  four  hydrophones  allow  the  computation  oi  the 
azimuth  and  elevation  angles  (spherical  coordinate  longitude  and  latitude)  of 
the  normal  to  the  wavefront  at  the  origin  of  the  local  coordinate  system. 
Then,  assuming  direct  path  propagation,  one  can  ray  trace  using  Snell's  law, 
[1],  starting  with  the  aforementioned  elevation  angle  and  utilizing  a  velocity- 
versus-depth  profile  for  the  speed  of  sound  in  the  water.  Finally,  the  time 
differential  between  the  source  pulse  at  the  target  vehicle  and  its  arrival  at  the 
array  is  used  to  stop  the  ray-tracing  algorithm  and  determine  the  location  of 
the  target  relative  to  the  array.  The  local  track  is  the  sequential  set  of  these 
estimated  positions. 

Each  array  in  the  system  operates  over  a  limited  radius.  As  the  target 
sojourns  through  the  range,  it  is  tracked  by  a  number  of  these  arrays.  See 
Figure  2  for  a  plan  view  of  the  Nanoose  range.  (The  zero  level  in  the  vertical 
is  taken  as  mean  sea  level.)  The  overall  path  is  constructed  by  transforming 
each  piece  of  local  track  to  the  coordinates  of  the  range  based  upon  the 


assumed  location  and  orientation  of  the  various  local  coordinate  systems. 
Discontinuities,  or  mismatches  occur  because  the  track  produced  by  one  array 
does  not  mesh  well  with  that  produced  by  a  neighboring  array  in  the  overlap) 
regions.     Such  mismatches  can  be  seen  in  Appendix  B. 

It  appears  that  there  are  several  sources  of  systematic  error  in  the 
operation  of  this  system.  The  question  of  individual  array  location  and 
orientation  has  been  treated  previously  [2].  The  present  report  addresses  the 
timing  synchronization  problem,  from  a  statistical  point  of  view.  That  is,  the 
pulses  received  at  the  arrays  are  timed  to  the  range  clock  with  great  precision. 
The  timing  information  must  be  transmitted  to  the  pinger  prior  to  the  release 
of  the  target  vehicle.  Although  there  is  no  engineering  reason  to  suspect 
noticeable  error  in  this  transfer,  some  data  exhibit  behavior  consonant  with 
such  an  interpretation.  Thus  we  build  a  mathematical  model  to  account  for 
such  a  source  of  systematic  error,  and  develop  statistical  methodology  for 
interpreting  the  data  in  the  light  of  the  model.  Indeed,  if  the  method 
provided  a  way  to  eliminate  mismatches  for  an  entire  run  for  a  single 
vehicle,  there  would  be  great  temptation  to  use  it  as  a  smoothing  filter. 

Generally  there  are  several  sequences  of  time  points,  called  point  count 
sets,  for  which  two  arrays  simultaneously  produce  track.  These  occur  for 
tracking  in  the  overlap  regions  indicated  in  Figure  2.  The  paired  tracking  data 
of  these  point  counts  are  called  "crossover  data."  It  is  the  crossover  data 
generated  as  a  result  of  the  target  vehicle's  entire  trip  that  provide  the 
evidence  suggesting  timing  synchronization  problems  and  the  data  base  for 
evaluating  the  use  of  such  a  model.  Each  version  of  track  in  a  crossover  data 
set  is  assumed  to  have  been  converted  to  the  (common)  range  coordinate 
system.  The  components  of  this  system  are  called  "downrange,  crossrange, 
vertical,"  or  sometimes  "centerline,  crossline,  vertical." 
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Figure  2.    Nanoose   Range 


The  model  is  developed  in  Section  2.  It  allows  for  both  a  timing  offset 
and  a  drift.  Estimation  methods  and  statistical  properties  are  developed. 
Section  3  contains  applications  to  real  data.  It  is  shown  that  the  offsets  are 
significant  and  the  drifts  are  negligible.  The  issue  of  the  reality  of  the  effects  is 
also  treated.  I.e.,  is  it  the  timing  offset  constant  for  all  crossover  regions  used 
in  a  given  run?  To  answer  this  question,  some  modeling  of  the  components 
of  variance  is  required  and  when  done,  a  statistical  analysis  is  performed.  It 
appears  that  the  true  cause  of  these  effects  is  something  other  than  a  timing 
synchronization  bias.  The  analyses  are  supported  by  some  characteristics  of 
the  noise  process,  and  these  are  presented  in  Section  4.  Conclusions  are 
summarized  in  Section  5.  A  number  of  appendices  are  included  that  contain 
supporting  data  and  information,  including  source  codes  for  computations. 

2.     MODEL  DEVELOPMENT  AND  STATISTICAL  PROPERTIES 

Figure  3  contains  a  mockup  illustrating  conditions  that  support  the 
consideration  of  a  timing  synchronization  correction.  It  shows  a  plan  view 
which  includes  the  radial  lines  from  the  arrays  to  the  estimated  track  points 
in  adjacent  overlap  regions.  Between  these  regions,  much  data  is  supplied 
only  by  array  A2.  The  analyst  does  not  see  the  radial  lines  on  his  screen,  nor 
does  he  see  the  black  dots.  He  sees  only  the  X's  and  O's  (no  distinction 
between  the  two)  and  no  mismatch  is  apparent.  (Mismatches  would  be 
apparent  however  for  track  pointing  in  a  different  direction.)  But  when  one 
pairs  up  the  radial  lines  by  common  time  points,  then  one  sees  that  the  two 
versions  of  track  lack  coherence  and  can  be  improved  by  stretching  the 
estimated  points  to  the  black  dot  positions.  This  can  be  achieved  by  a  single 
constant  adjustment  to  the  transit  time  values  in  the  ray  tracing  algorithm. 
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The  figure  also  helps  one  to  imagine  the  distinction  between  the  effects  of 
an  array  location  error  model  and  a  timing  adjustment  model.  If  the  former 
were  applied  to  the  situation  in  Figure  3,  at  least  two  arrays  would  have  to  be 
moved.  But  any  such  decisions  cannot  be  made  in  isolation.  The  ultimate 
goal  is  the  simultaneous  improvement  of  coherence  in  all  overlap  areas  and, 
of  course,  discontinuities  must  not  be  created  at  other  places. 

A  crossover  data  set  is  a  set  of  matched  pairs  of  three  vectors 


X(t)  = 


X^t) 
X2(t) 

lx3(t)j 


Y(t)  = 


[Yi(0~ 
Y2(t) 
Y3(t)J 


t«l,...,T 


(1) 


representing  two  versions  of  the  same  vehicle  track  for  a  common  set  of  time 
values  (point  counts),  T  in  number.  The  array  that  produced  X(t)  is  located  at 
a  and  the  one  that  produced  Y(t)  is  located  at  (3.  Thus  these  two  versions  of 
track  can  be  represented  in  their  local  coordinate  systems  as 

%(t)  =  X(t)-oc  ri(t)  =  Y(t)-P  (2) 

Timing  synchronization  error  corrections  may  be  viewed  as  either  stretching 
or  contracting  the  vectors  £,(t)  and  rj(t)  by  the  same  (time)  amount.  The  effect 
of  such  corrections  is  not  constant,  but  depends  upon  (i)  the  speed  of  sound  in 
water  at  the  depths  X3(t)  and  Y3(t);  (ii)  the  elevation  angles  of  the  ray  traces  at 
these  times  and  depths. 

The  magnitudes  of  the  time  adjustments  are  small  and  the  effects  can  be 
described  using  first  order  terms 


(3) 


£(t,g(t))  =  5(t)  +  g(t)a(t) 

T|(t,g(t))  =  T](t)  +  g(t)b(t) 

where  g(t)  is  a  scalar  function  of  time  representing  the  total  adjustment  for 
timing  offset  (5)  and  drift  (m) 


g(t)  =  5  +  m(t-to) 


(4) 


for   a    conveniently   chosen   central   time,   to;  and  a(t),  b(t)  are  the  scaled 
directions  of  stretch 


rcos(6(t))cos«)(t))>' 


a(t)  =  v(t)  i  cos(9(t))sin(0(t)) 
sin(9(t» 


(5) 


i.e.,  v(t)  is  the  speed  of  sound  at  depth  ^(t),  6(0  is  the  elevation  angle  of  the 
ray  trace  from  a  at  ^(t),  and  0(t)  is  the  azimuth  of  ^(t)  from  a.  The  vector  b(t) 
is  defined  similarly  for  the  track  r|(t)  and  its  origin  p. 

It  is  easily  seen  that,  for  given  5,  m  and  to,  the  corrected  versions  of  track 
in  the  (common)  range  coordinate  system  are 

X(t,g(t))  =  X(t)  +  g(t)a(t) 
Y(t,g(t))  =  Y(t)  +  g(t)b(t) 

The  stat:stical  estimation  problem  is  to  chose  8  and  m  so  that  the  two 
versions  of  track  agree  as  well  as  possible.  The  least  squares  approach  is 
adopted.    Using  the  notation 
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Ave[W(ff  =i]TW'(0W(0 

1  '  f=i 


(?) 


we  set  up  the  objective  function 


Q  =  Aye\\x{t,g(t))-Y{t,g(t)f 

=  Ave|X(f)-Y(/)||2  +  Ave[g(t)f\\a(t)-b(tf  (8) 

+2Ave  g(tla(t)-b(t)][X(t)-Y(tj\ 

and  find  the  values  of  5  and  m  that  minimize  Q.  When  the  two  partial 
derivatives  are  set  equal  to  zero  it  is  seen  that  there  is  convenience  in 
choosing  to  so  that 

Ave(t-t0)|K0-Wl2  =  0  (9) 

That  is 

tQ  =  Avet\\a(t)-b(tf/Ave\\a(t)-b(tf 

t  t 

This  done,  the  normal  equations  may  be  expressed  as 

8Ave\\a{t)-b{tf  =  -Ave(a{t)-b{t))[X{t)-Y(t)] 


(10) 


mAve(t-t0)%(t)-b(t)f=-Ave{t-t0)[a(t)-b(t)][X(t)-Y(t)] 


and  one  can  readily  solve  explicitly  for  the  minimizing  values  <5  and  m. 
Positive  values  of  8  are  associated  with  stretching,  negative  with  contracting. 
If  we  let  Q  be  the  minimizing  value  of  Q,  it  is  useful  to  establish  the 
decomposition 

Q  =  Q  +  (S  -  sf  Ave|fl(r )  -  b(tf  +  [m  -  mf  Ave(f  - 10 )\\a{t)-  b(tf        (11) 


10 


Before  justifying  (11)  it  is  convenient  to  shorten  the  writing:  let  \\7(t,6(t)) 
=  X(t,g{t))-Y(t,g(t));  c{t)  =  a{t)-b{t);g{t)=8  +  m(t-tQ).  Then  the  statement 
says 

Aye\w(t,g(t)f 

=  Aje\w(t,g(t)f  +  [S  -  sf  Ave||c (Ol2  +  (m  -  m)2  Ave(f  -  f0)||c(0|2  (12) 

and  the  result  is  justified  by  use  of  the  three  orthogonality  relationships 


Ave 

t 


W{t,0)+  5c(t)  +  m(t-t0)c(t)   c{t)  =  0 


Ave 

t 


W(t,0)  +  5c(t)  +  m{t-t0)c(t)]{t-t0)c(t)  =  0 
Ave{t-t0)c'(t)c(t)  =  0 


(13) 


which  in  turn  are  established  using  the  normal  equations  (10),  and  (9). 

The  significance  of  the  offset  and  drift  parameters  can  be  judged  if  we 
develop  the  means  and  variances  of  8 ,  m.  Letting  E  denote  the  mathematical 
expectation  operator,  we  begin  with  the  assumption 

E[W(t,g(t)]  =  0  (14) 

which  embraces  the  idea  that  the  two  corrected  versions  of  track  produce 
common  track  without  systematic  error. 
Further,  let 


K.  = 

6 


^c'(t)c(t)\    ;        Km  = 


Z(t-t0)2c'(t)c(t) 


] 


(15) 


and  use  (13)  and  (14)  to  show  that  the  estimators  5  and  in  are  unbiased. 


]  I 


£(<S)  =  -K*5V(04W(''°)]  =  KS^V)[S  +  m(t  -  t0)]c(t)  =  8 
t  i 

E(rh)  =  -Km^(t  -  t0y(t)E[W(t,OJ\  =  Km^(t-  t0)c'(tp  +  m(t  -  t0j)f(t)  =  m 

i  t 

To  develop  variances,  we  assume  that  the  positive  lag  covariances  are 
zero;  support  for  this  appears  in  the  section  on  noise  characteristics.  Let  M 
represent  the  (zero  lag)  covariance  matrix.    (See  Appendix  C  for  estimates). 

M  =  E[w(t,g(t))W'(t,g(tj\  (16) 

Then  one  easily  represents 

var(5)  =  Kj£c'(t)Mc(0  (17) 

t 

wax{m)  =  K2m^{t-t0)2c\t)Mc\t)  (18) 

cov[8tm)=K5KmY,{t-tQy(t)Mc(t)  (19) 

l 

If  M  is  proportional  to  the  identity  matrix  then  this  last  term  is  zero;  the  other 
two  terms  simplify  immensely;  and  a  standard  regression  development  can 
be  used.  But  the  study  of  noise  characteristics  does  not  support  this.  On  the 
other  hand  the  vectors  {c(t)}  do  not  change  much  with  t  and  this  has  the 
tendency  to  render  (19)  to  be  small.  The  reason  for  this  stability  is  that 
crossover  data  occurs  only  at  the  greater  distances  from  the  sensing  arrays;  the 
azimuth  and  elevation  angles  and  the  layer  sound  speeds  do  not  change 
much. 
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It  appears  that  the  matrix  M  changes  with  the  crossover  data  set. 
Methodology  for  estimating  it  appears  in  Section  4.  Estimates  of  the  M 
matrices  appear  in  Appendix  C. 

3.     DATA  ANALYSIS 

The  data  consists  of  8,  m  ,  and  supporting  values  for  69  segments  of 
crossover  track  collected  over  four  separate  days  with  two  (temporally  serial, 
not  concurrent)  target  vehicles  (runs)  per  day.  The  information  is 
summarized  in  Table  1.  Missing  variance  estimates  indicate  either  a  data 
shortage  or  outlier  problems.  In  a  few  cases  the  tracks  were  curved,  and  the 
straight  line  model  is  not  adequate.  The  units  are  milliseconds  for  8,  and 
milliseconds  per  point  count  for  m.  One  millisecond  translates  to  about  4.5  to 
5  feet  of  distance.  The  ratios  of  means  to  standard  deviations  are  used  to 
judge  whether  the  effects  are  significantly  different  from  zero.  Virtually  all  of 
the  offsets  and  some  of  the  drifts  are  significant  although  the  latter  are  not 
strongly  so.  The  "r"  column  contains  the  correlations  between  8  and  m  . 
They  are  insignificant.  A  further  search  for  large  scale  drift  was  made  by 
plotting  8  against  the  crossover  central  time  to  for  each  run.  They  appear  in 
Figure  4.  If  a  smooth  signal  were  discernible  then  we  would  have  a  way  to 
connect  the  5  values  that  appear  in  each  column  of  Table  2.  But  they  are 
chaotic  and  provide  no  incentive  to  continue  any  concern  about  drift.  It  is 
concluded  that  the  offsets  are  significant  and  the  drifts  are  not.  The  latter  will 
be  dropped  from  further  consideration.  Some  graphical  examples  of  the  effect 
our  timing  corrections  have  been  selected  and  appear  in  Appendix  B. 
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Table  1:    Offsets,  Drifts,  Signal  to  Noise  Ratios 


Date 

3/2  3/89 


Vehicle 

MK  30       1 

1 
i 


3/2  3/8  9        MK  30 
5/1  0/89        MK  27 


5/10/89        MK  27 


6  /  6  ,'  8  9 


MK  27 


6    6    89 


MK  27 


7/21  /88        MK   27 


7/21/88         MK   30 


37 
74 
69 
78 

1  5 

10 
24 

54 
52 
04 
44 

.24 
43 
86 
24 
06 
99 
27 
89 
49 
66 
90 
93 
09 
65 

.87 
37 
38 

.36 
18 
'6 
91 
1  5 
68 
31 
95 
72 

0  3 

1  4 
9  1 
6  9 
1  1 
0  9 
9  0 
4<r 
99 

0  4 
3  2 
07 
^5 
24 
9? 
47 
16 
29 
45 
3" 

1  5 
96 
10 
6  9 
C  2 

80 

42 
2~ 
17 
£  4 
61 

4~ 


0461 
0431 
0925 
0464 
0504 
0802 
0480 
0479 
0584 
0283 
0192 
.0401 
.0270 
0266 
0228 
0197 
.0104 
0180 
0868 
0194 
0270 
0189 
0258 
0406 
0243 


0131 
0209 
0181 
0149 
0096 
0281 
0336 
01  16 
0269 
0125 
0136 
0376 
0251 
0151 
0232 
0434 
0195 
0233 
0156 
0202 
0181 
0217 
0300 
01  15 
0147 

0128 
0144 
01  12 
0245 
0308 
0192 
0171 
0189 
0099 
0209 
0337 
0176 
0251 
0496 
0259 
02<~ 


fro, 


'/• 

29  82 
40.28 
18  29 

16  87 

2.97 

1  23 

25  96 

73  85 

43  1  7 

107.40 

178  88 

80  68 

90  18 

69  80 

54  56 

104  59 

190  86 

126  16 

21  72 

180  1  1 

135  56 

259  57 

1  13  27 

51  50 

117.13 


257  71 

64  92 

120.55 

146  48 

510  80 

1  1  1  97 

79  78 

371  65 

184  19 

21  7  40 

297  47 

30  32 

76  25 

178  04 

91  30 

7115 

148  86 

106  78 

127.36 

51  69 

72  74 

95.45 

91  58 

195  77 

129  27 

90  47 
19  86 
217  86 
55  85 
69  77 
50  18 
63  97 
36  48 

2  07 
38  31 
12  35 
15  46 

6  69 
10  88 
23  51 
1  5  94 


m 

00168 
00831 
01  559 
00245 
.01304 
00150 
001  1  9 
00395 
00224 
00993 
00250 
00007 
00187 
00029 
00408 
00437 
00067 
.00196 
00001 
01062 
00044 
00443 
00545 
00038 
00094 
.00050 
.01074 
00257 
.00276 
.00195 
00192 
00022 
00202 
01680 
00189 
.00144 
00343 
00242 
00197 
.00401 
001  43 
00091 
00596 
00231 
.0021 1 
00192 
00075 
.00003 
00060 
.00020 
00160 
00244 
.00746 
00287 
.00232 
00069 
01  149 
.00069 
00019 
.00206 
00233 
.00190 
.00147 
00355 
00035 
00024 
001  71 
00509 
"'  ■  23 


001627 
002217 

.007024 
001 588 
001438 
005509 
001645 
0021 19 
002136 
001 142 
000866 

.002747 

.001375 
000825 
000630 
000662 
000364 

.000770 
000567 

.000746 
001268 

.000689 

.001358 
001887 
000845 


000458 
000662 
000628 

.000479 
000337 
001001 
003073 

.000377 

.001 136 
000434 

.000464 
001555 
001 141 

.000486 
000783 
003137 
000677 
000692 
000541 
000681 

.000631 

.000725 
001845 
000398 
000510 

001310 
000501 
000392 
000853 
.002324 
000662 
000595 
000656 
000346 
000730 
001 173 
00061 1 
000952 
003582 
000897 
001380 


m  /C„ 


1  03 

3  75 

2  22 

1  54 

9  07 

27 

72 

1  86 
1.05 
8  70 

2  89 
02 
36 

.35 
48 
60 

1  84 

2  54 
01 

14  23 

35 

6  43 

4  02 
.20 

1.11 


5  62 
4  17 

3  1  1 

4  02 
64 
02 
47 
02 
27 


9' 
22 
27 

5  2 
93 

1  17 
1  90 
42 
06 
55 
1  0 

0  5 

6  3 

1  1 
01 
79 


2  19 
4  64 

1  76 

13  46 

.30 

29 


46 
6  4 
49 

01 
02 
57 

2  6 
48 
67 
34 


00 
00 
06 
02 
02 
01 
01 
06 
1  1 
02 
06 

.05 
00 

.00 

.05 
02 

.02 

.03 
01 
03 

.00 
0  1 
00 
01 

.07 


00 

06 

07 

.03 

.02 

.05 

.02 

08 

.04 

.02 

.02 

.01 

.07 

05 

03 

.02 

03 

.04 

02 

.04 

01 

02 

00 

.00 

.03 

01 
01 
00 
02 
02 
03 

.02 
02 
06 
00 
03 
07 
03 

.03 
02 
04 


to 

405  13 

422.72 
4171  09 

439  71 
2012  80 
3080  16 
2531  21 

422  64 
1605  22 
3107  58 
3120  22 
6936  70 
3116  86 
3781  03 
5471  41 
4720.56 
8060  02 
5992  35 
961 5  69 
1733  39 
1751  02 
0196  44 
1  747.85 
2366  96 
2980  01 
2354.91 
5316  80 
2215  15 
3585.00 
4478  00 
5873  92 
1235  86 
6480  62 
1952  64 
0562  53 
441  1  90 
2124  12 
2709  75 
4979  22 
6197  42 
2044  01 

357  26 

421  96 
4135  53 

487.61 
1781  82 
2505  37 
1991  06 

352  15 

377  91 
4440  01 

469  17 
1062  61 
1772  98 
2068  58 
3137  66 
1 108  15 
1 152  98 
6658  19 

217  45 
4360  14 
1498  22 
4667  48 
0430  20 
3793  99 
5141  88 
2031  94 
301  1  04 
5697  10 


4  2 

3  3 
1  9 
45 

5  7 
22 
49 
37 

4  7 

4  3 

3  6 
25 
34 

6  6 
59 
50 
50 
40 
50 
45 
37 
45 
33 
37 

5  0 
50 
1  5 
50 
55 
50 
54 
5C 
46 
1  9 
50 

4  • 

50 
50 
40 

:•  5 
51 
50 
24 
50 
50 
50 

5: 

50 

5  2 
28 
50 
50 
1  5 

1  7 

5: 

50 
5  0 

2  3 
5  0 
50 
50 
50 
50 
50 
50 
27 
24 
50 
37 
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Table  2:    Estimates  of  the  Timing  Offset 


4.91    (i)    4-95   (k> 
4.89  (k) 


4.31    (k) 


4.03  (e) 


3.0- 


2.5- 


2.0  — 


1.5- 


1.0- 


0.5  — 


0   — 


-0.5  — 


-1.0 


3.86  (e) 
3.66  (b) 


3.53   (b) 


3.44   (b) 
3.23   (k) 
3.04  (a) 

3.38  (h) 
3.14  (g) 
2.92   (C) 

2.85  (e) 

2.52  (e) 

2.43  (c) 

2.27  (g) 

2.37  (h) 
2.17    (f) 

2.05  (h) 

2.17    (f) 

1.98    (j) 

2.08  (d) 

1.88  (g) 

1.73  (b)                    185  (d) 

1.69  (k) 

308  (b) 
2.89   (k) 


2.71  (c) 
2.69  (f) 
2.68  (a) 


2.74   (b) 


1.37  (a) 


0.78   (C) 


0.09  (g) 

•0.14    (h) 


■1.24    (f) 


' 1 

3   2'    89 


1.91    (f) 


2.48   (C)  2.44 
2.24 

2.11   (a)2Q7 
1.99  (e)  ^go 


1.24    (f) 


1.35    (f) 


1.14    (f) 


1.31    (h) 


1.04    (f) 


1.47 


1.16 


(k) 
(a) 

(c) 

(d) 
(e) 


2.14    (b) 


0  28  (h) 


1.36  (a) 

1.09  (c) 
0.96  (k) 

0.68  (e) 
0.60    (j) 

0.46  (g) 

0.27    m 
0.16    (f) 

-0.02    (e) 

-0  41    (h) 
-0.53    (f) 

-0.80    (h) 


j 


'89      U 


6/8  9 


7/21    88 


3/23/89  -MK  30 


3/23/89 -MK  30 


3  6  ._ 

• 

3  4  • 

3.2  • 

1  . 

2-1  • 

2.6  ■ 


0.6  08  1  1.2  1.4 

t.xlO"3 


5/10/89 -MK  27 


5/10/89  -MK  27 


1 


5  6  7  6 

t„xl0'3 


6/6/89 -MK  27 


6/6/89 -MK  27 


0  0.5  1 


t„xl0 


t.xlO" 


721/88-MK27 


7/21  /88  -  MK  30 


Figure  4.   Offset  vs.  Central  Point  Count  for  the  eight  runs 


Since  the  <5  values  are  significant,  it  is  important  to  discover  what  they 
really  represent.  They  have  been  labeled  timing  synchronization  offset  errors 
and,  if  that  labeling  is  a  valid  physical  description,  then  each  run  (vehicle)  on 
each  day  should  have  its  own  unique  value  and  this  common  value  could  be 
used  to  correct  all  track,  not  just  crossover  track.  What  follows  is  an  analysis 
of  this  point. 

Table  2  contains  a  graphical  positioning  of  the  sixty-nine  S  values 
according  to  the  eight  runs.  There  were  eleven  array  pairs  involved  in  this 
process  and  they  are  marked  with  the  letters  a,  b,  ...,  k.  (The  correspondence 
with  the  overlap  regions  marked  in  Figure  2  is  given  in  Table  3.)  Study  of 
Table  2  shows  considerable  "within  run"  variability,  and  one  is  inclined  to 
doubt  the  reality  of  the  constant  offset  interpretation. 

TABLE  3.  IDENTIFICATION  OF  ARRAY  PAIRS  TO  RECONCILE  TABLE  2 

WITH  FIGURE  2 


a      (1,11) 

e      (3,4) 

i       (5,56) 

b      (1,2) 

f       (43) 

j       (5,14) 

c      (2,11) 

g      (5,6) 

k      (15,16) 

d      (2,3) 

h      (4^5) 

In  passing  we  also  note  that  the  S  values  for  each  array  pair  appear  to 
have  some  temporal  coherence,  and  hence  the  assignable  causes  may  be 
related  to  this,  but  that  is  an  issue  for  another  time. 

We  proceed  to  model  the  components  of  random  error  variance  and 
develop  a  test  statistic  for  examining  whether  or  not  all  5  for  the  same  run 
can  be  viewed  as  constant.  The  {<5,}  are  modeled  as  being  affected  bv  the  dav 
(water  depth  velocity  profiles  change  from  day  to  day),  the  run  (the  second 
run  is  a  different  vehicle  operating  later  in  the  day),  and  the  array  pairs 
generating  the  crossover  data.    The  array  pairs  are  assumed  to  produce  fixed 
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effects,  with  values  a,  for  j=l,  ...,  11;  and  the  day  and  run  effects  are  assumed 

2  2 

to  be  random  with  zero  means  and  variances  O]  and  o2  respectively.    There  is 

2 
also  an  experimental  (residual)  error  having  zero  mean  and  variance  a  . 

The  result  is  the  mixed  model 

8  =  Xa  +  Zp  +  e  (20) 

where  5  is  the  69  component  column  vector  of  6  estimates;  X  is  a  69  x  11 
matrix  of  zeros  and  ones  (incidence  matrix)  relating  the  8  component  to  the 
array  pair;  a  is  the  11  component  vector  of  array  pair  fixed  effects;  Z  is  a  69 
rowed  partitioned  incidence  matrix 

Z  =  (Z1/Z2)  (21) 

with  Z]  having  4  columns  relating  the  8  component  to  its  day  and  Z2  having 
8  columns  relating  the  8  component  to  its  run.     The  vector  P'  =  (p.,  p2) 

represents  the  random  variables 

pj  =  (Pn,  ...,  (314)       independent      N(0,o1) 

2 
P2  =  (P21/  ••■/  P28)      independent      N(0,O2) 

and  the  residuals  are 

e'  =  (ei,  ...,  e69)        independent      X(0,a2) 

Further  the  vectors  pi,  P2,  e  are  assumed  to  be  mutually  independent. 

We  will  be  applying  the  model  (20)  separately  to  each  of  the  eight  runs  for 

purposes  of  estimating  the  fixed  effects  (a)  for  each  of  those  runs.     Such 

2       2 

estimates  will  require  values  for  the  variance  components  (o2,  c1#  02).  We 
prefer  to  estimate  them  but  once,  by  pooling  all  oi  the  data.     This  can  be 


is 


accomplished  separately  from  the  estimation  of  fixed  effects  by  using  the 
restricted  maximum  likelihood  method:  the  unequal  block  size  version  has 
been  treated  by  Patterson  and  Thompson,  Ref.  [4]  and  is  used  here.  (See 
Appendix  A  for  a  description  of  the  algorithm.)   The  results  are 

61=  0.3571         02  =  0.7207       a2  =  0.2430  (22) 

and  the  units  are  seconds  squared.  It  is  noted  that  the  run-to-run  variance  is 
double  that  of  the  day-to-day  variance,  which  in  turn  is  larger  than  the  error 
variance  by  a  ratio  of  about  seven  to  five.  But  the  degrees  of  freedom  for  the 
first  two  variances  are  small  and  the  estimates  may  not  be  very  reliable. 

Using  the  model  (20)  one  sees  that  the  covariance  matrix  of  the 
observables  is 

H  =  o2In  +  ZrZ'  (23) 

where  In  is  the  identity  matrix  of  order  n,  V  is  diagonal,  the  first  four  values 

2  '  2 

being  a-j   and  the  next  eight  02;  and  Z  is  partitioned  (Zi,  Z2)  as  before.    For  the 

full  data  set  n  =  69,  but  recall  that  our  goal  is  to  check  whether  the  fixed  effects 
(offset  estimates)  can  be  viewed  as  constant  within  each  run.  The  plan  is  to 
estimate  the  fixed  effects  for  each  run  and  test  whether  they  can  be  viewed  as 
constant  (over  the  array  pairs). 

To  do  this  we  proceed  as  follows.  The  5  values  for  the  kth  run  can  be 
identified  by  the  ones  in  the  kth  column  of  Zo,  k  =  1,  ...,  8.  The  array  pairs 
involved  (not  necessarily  all  eleven)  are  identified  by  ones  in  the  X  matrix 
restricted  to  the  kth  run.  The  Z  matrix  is  restricted  in  a  like  fashion  and  the 
covariance  matrix  (of  order  n,  the  number  of  crossover  data  sets  in  the  kth 
run)  has  the  same  structure  as  (23)  with  the  reinterpretation  oi  inputs.  Of 
course  the  estimates  (22)  must  also  be  input. 


19 


Let  K  be  the  number  of  array  pairs  involved  in  the  kth  run.     Use  the 
Aitken  estimator  [3] 

ak  =  (xkH?Xk)~lXkH?8k  (24) 

and  its  covariance  matrix 

vai(ak)  =  [xkHk%]~1  =  V~^  (25) 

and  the  subscript  k  modifies  the  previous  definitions  of  symbols  so  that  only 
the  kth  run  is  involved.  If  the  null  hypothesis  is  true  (i.e.,  the  ak  all  fall  on 
the  main  diagonal  of  the  K  dimensional  space),  then  the  maximum 
likelihood  estimate  of  the  common  value  is 

« =ii^M/)/ii^  (26) 

It  also  follows  from  the  normality  assumptions  that  the  quadratic 

(ak-a)V(ak-a)      ~      Chi  Square  (K-l)  (27) 

and  this  statistic  is  the  weighted  distance  of  the  components  of  &  from  the 
main  diagonal  (i.e.,  the  constant  fixed  effect  that  is  the  same  for  all  array 
pairs).  Thus  the  null  hypothesis  should  be  rejected  when  this  distance  is  too 
large. 

The  numerical  results  for  the  eight  runs  are  contained  in  Table  4.  The 
column  marked  p*  contains  the  p  values  (empirical  significance  level) 
corresponding  to  the  test  statistics  listed  under  the  "distance"  column.  Seven 
of  the  eight  values  indicate  rather  rare  events  and  the  eighth,  p*  =  0.14,  is 
associated  with  one  degree  oi  freedom.     Low  degrees  of  freedom  tend  to 
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dampen  the  chances  of  finding  significant  results.    The  final  column  is  the 
ratio  of  the  distance  (27)  to  its  standard  deviation  under  the  null  hypothesis. 

We  assert  that  the  array  pair  effects  are  not  constant  and  that  there  are 
other  sources  of  systematic  error  that  dominate  this  process. 

TABLE  4.  TESTING  THE  EIGHT  RUNS  FOR  CONSTANT  FIXED  EFFECTS 


a 

K-1 

Distance 

P* 

distance/SD 

0.61 

6 

30.17 

3.6x1 0-5 

8.71 

3.03 

1 

2.12 

0.14 

1.51 

2.35 

8 

17.71 

0.06 

4.43 

3.09 

9 

51.67 

5.2x1a-8 

12.18 

3.05 

4 

41.48 

2.1x10-8 

14.67 

2.13 

6 

14.16 

0.03 

4.09 

1.79 

7 

18.13 

0.01 

4.85 

0.46 

8 

30.44 

1.8xl(H 

7.61 

4.     CHARACTERISTICS  OF  NOISE 

In  order  to  study  the  stochastic  nature  of  the  tracking  errors  we  fit  straight 
lines  to  the  track.  (This  was  done  even  for  track  that  did  not  appear  linear  via 
a  visual  scan  of  graphical  output.  The  exceptional  cases  are  marked  in 
Appendix  C.)  It  is  assumed  that  the  target  vehicles  speed  is  constant.  The 
result  is  a  set  of  deterministically  spaced  points  that  fall  on  a  straight  line  in 
three  space.  From  these  we  can  produce  residuals  and,  assuming  local 
stationarity,  study  the  noise  structure. 
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The  basic  (unadjusted  for  constant  speed)  straight  line  is  found  by 
principal  components.  Let  X(t)  be  used  for  X(t,0),  eq.  (6)  and 

5f-l£X(0  (28) 

'    t 

We  seek  a  projection  (direction)  p  =  (pi,  p^  P3)' 

Z(t)  =  j^plXl(t)  =  p'X(t)  (29) 

l 

that  will  maximize  the  variance  of  (Z(t)} 

°Z=^l[Z(t)-zf  =  P'Cxp  (30) 

1    t 

where 

Cx=^(X(t)-X)(X(t)-X)  (31) 

t 

is  the  covariance  matrix  of  the  track  data,  X(t). 

Since  the  vector  p  is  a  set  of  direction  numbers,  we  adopt  the  usual 
constraint,  2*1  P;  =  *•  (This  will  keep  oz  finite.)  It  is  easily  shown  that  the 
three  solutions  to  the  eigen  problem 

Cxp  =  Xp 

provide  us  with  an  orthonormal  basis 

P  =  {PuP2fP3]  (32) 

2 
where  P,  is  the  eigen  vector  corresponding  to  X\  and  X\  >  X2  >  X3.   Since  o"z  = 

2 
p'Cxp  =  ^p'p  =  X,  and  our  goal  is  to  maximize  G^,  we  choose  the  first  eigen 

vector  for  use. 
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The  set  of  values 

Z1(f)  =  P1'(X(0-X)  (33) 

represents  the  succession  of  values  after  the  vectors  {x(f)-X}  are  projected 

onto  the  line  of  the  first  principal  component.  These  projections  are 
orthogonal  projections  and  may  not  well  represent  the  relative  positions  of  a 
target  vehicle  heading  in  the  direction  Pj  at  constant  speed.  Accordingly  it  is 
appropriate  to  make  an  adjustment  in  the  (Z](t)}  to  account  for  constant 
speed.   The  collection  of  point  counts, 

tl<t2<...<tTl  (34) 

all  differ  by  integral  multiples  of  some  base  value,  A,  representing  distance 
traversed  per  point  count.  Let  us  perform  a  simple  linear  regression  of  the 
Z](t)  on  the  times  (34).   The  fitted  values  are 

Z(t)  =  a  +  b{t-t)  (35) 

where  a  =  -^iZl{t)  and  b  =  £  (f-fjZ^f)  /£  (f  -  f ).   Now    the    successive 

values  of  Z(t)  differ  by  integral  multiples  of  a  converted  base  value,  bA,  and 
represent  the  constant  speed  progression  of  the  vehicle  in  the  direction  P]. 

It  remains   to  represent  the  distances   Z(t)   in   the  original   coordinate 
svstem.  Let 

R(t)  =  [z(t),0,0)'  (36) 

be  the  estimated  position  of  the  vehicle  in  the  basis  P,  and  then 

Xd{t)  =  X  +  PR(t)  (37) 
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will  be  the  straight  line  in  the  original  coordinates.    Similarly,  when  (Y(t)}  is 
put  into  this  algorithm,  we  produce  Ys](t). 

Now  we  are  positioned  to  estimate  the  covariance  matrices  of  the  noise 
processes 

Dz=7X[*)-xs,W][x(')-xs,W]' 

t 

oy  =  ^[y(t)-y5i(t)}[y(t)-Ysl(t)]  (38) 


y    T 


t 


Dxy=^[X(t)-Xsl(t))[Y(t)-Tsl(t)] 
T   t 

Computational  work  shows  considerable  variability  in  these  matrices  as  one 
changes  the  day,  the  run,  and  the  array  pair.  Also  the  cross  correlations  are 
mostly  different  from  zero.   The  matrix  M,  eq.  (16),  is  estimated  by 

M  =  DX  +  Dy-Dxy-Dxy  (39) 

and  used  in  (17),  (18),  and  (19).  The  quantities  (38)  and  (39)  appear  in 
Appendix  C,  and  illustrate  their  variable  nature. 

Another  immediate  use  of  the  noise  processes  is  to  look  at  the 
autocorrelations.  In  addition  to  the  three  individual  components  we  define  a 
"noise  displacement"  process 

l 


d{t)=h(x,(t)-x5lt(t)f 


(40) 


(Of  course  the  same  is  done  for  (Y(t)}.) 
The  autocovariances 
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T—h 

R(h)  =  ^t(d(t)-d)(d(t  +  h)-d)  (41) 

for  (d(t)}  and  for  the  individual  components  {Xj(t)-Xsi(t)}  are  computed  and 

normalized.    Plots  oi  R(h)/R(0)  vs  h  appear  in  Appendix  D.    It  is  typical  that 
they  become  and  remain  in  a  general  level  of  "static"  for  h  >  1. 

5.     SUMMARY 

The  paper  contains  methodology  for  estimating  the  presence  of  timing 
synchronization  error  and  judging  its  significance  in  the  framework  of  our 
short  baseline  array  underwater  test  range.  The  methodology  was  applied  to 
real  data.  Some  plots  illustrating  the  effects  appear  in  Appendix  B.  One  must 
view  them  with  care  because  the  coordinate  scales  are  so  different.  It  is  rather 
tvpical  that  some  unresolved  systematic  error  is  exposed  in  the  side  view. 

Analysis  of  the  results  does  not  support  the  timing  synchronization  error 
model  as  accounting  for  the  discrepancies.  Other  sources  of  systematic  error 
must  be  unmasked  first.  Possibilities  include  the  orientation  of  the  arrays  [2], 
biases  in  the  raytracing  inputs  [3],  and  perhaps  a  temporal  or  spatial  gradient 
in  the  depth  velocity  profile. 

Appendices  C  and  D  contain  information  about  the  second  order 
properties  of  the  three  dimensional  noise  process  and  about  time  lag 
correlations. 
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APPENDIX  A.  PATTERSON  AND  THOMPSON  ALGORITHM 

The  purpose  of  this  appendix  is  to  document  the  computational  method 
for  estimating  variance  components  that  is  but  implicitly  described  in  the 
Patterson  and  Thompson  paper  [4].  We  use  the  notation  of  that  paper. 
Specifically 

y  =  Xa  +  e  (A.l) 

where  the  observeable  y  is  an  n  component  column  vector,  X  is  an  n  by  t 
matrix  of  rank  t  determined  by  the  allocation  of  treatments  to  units,  a  is  a  t 
vector  oi  fixed  effects,  and  e  is  a  mean  zero  n  vector  of  normal  random 
variables  with  covariance  matrix 

V  =  o2H       H  =  ZTZ'  +  I.  (A.2) 

The  matrix  I  is  the  identity  of  order  n,  Z  is  the  n  by  b  design  matrix  for  c  block 
factors,  and  I"  is  a  diagonal  matrix  containing  the  variance  components 
relative  to  the  basic  error  variance  o2.  I.e., 

c 
H  =  I  +  XZpZp'Yp  (A.  3) 

P=l 

and  each  block  design  matrix,  Zp,  is  n  by  bp;  each  plot  having  exactly  one  level 
in  each  Zp,  p  =  1,  ...,  c,  the  variance  of  the  pth  block  is 

2  ■> 

op  =  Yp<^       for  p  =  1,2,  ...,c  (A.4) 

and  Z  in  (A.2)  has  the  partitioned  form 

Z  =  (Z1,Z2/...,ZC).  (A.5) 
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The  diagonal  elements  of  V  are  b>i  consecutive  yi's,  hi  consecutive  72%  •••,  and 

c 

bc  consecutive  7c's,  with  X^p  =  b. 

P=i 

Out  goal  is  the  estimation  of  71,  72,  ...,  yc,  a2  by  the  restricted  maximum 

likelihood  method  for  the  unbalanced  case.   The  steps  follow.   Let 

S  =  I-X(X'X)"1X'  (A.6) 

W  =  Z'SZ+r"1  (A.7) 

The  algorithm  is  an  iterative  one,  and  initial  values  for  the  71,  ...,  7C  are 
needed  to  develop  W  in  (A.7).  The  b  bv  b  matrix  VV  and  its  inverse  can  be 
partitioned  according  to  b],  D2,  ...,  bc.  So  doing  allows  both 

(3  =  W-1Z'Sy  (A. 8) 

u  =  r^-r'V'Y"1  (A.9) 

to  be  viewed  as  partitioned.   I.e., 


(  '     '  '  \ 

I  ) 

U  =  {Ujj)    i,  ]  =  !,... ,c 


(A. 10) 


and  Uij  is  a  bj  bv  b,  matrix. 

Now  we  are  positioned  to  define 

R  =  y'Sy-y'SZp 

Pp  =  PpPP/Yp  (A11) 

and  Ep  =  trace  (Upp)  for  p  =  1,  •-.,  c.    Then  define  the  modified  information 
matrix  (fj)},  of  order  c+1, 
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fij  =  trace  (UjjUji)      for  i,  j  =  1,  ...,  c 
fp,c+i  =fc+hp  =  Ep       forp  =  l,...,  c 

fc+l,c+l  =  n~t 

and  let  {f'i}  be  its  inverse. 

Finally,  the  estimator  update  equations  are,  using  k  =  c+1 


G2  =  Y.f^B^fkkR 

y,  =  y    ;  +  — 


(A.12) 


7i  =  7, 


o~ 


(A.13) 


for  i  =  l,...,c 


and  the  new  {/,-}  can  be  placed  in  (A. 7)  to  start  the  next  iteration.  No  initial 
value  for  o2  is  required,  but  it  must  be  updated  at  each  iteration.  The 
algorithm  should  be  stopped  when  (A.13)  is  stable.  The  variance  components 
are  computed  from  (A. 4) 

The  application  in  the  present  report  has  c=2,  with  b]  =  4  and  b2  =  8.   It  is 
helpful  to  record  the  partitioning  and  inversion  of  W,  (A. 7). 


W  = 


/-2J^\      Z.20Z.2 

wn    w12 

w21    w22 


+ 


0  I2/y2 

symmetric. 


(A.14) 


W"1  = 


w11    w12' 
w21   w22 


can  be  computed  from 


(A.15) 
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w22=[w22-vv21Wn1vv12] 
w12  =  -wrfwnw22 
w11  =  [/1-w12w21]w1-11 


1 


(A.16) 
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APPENDIX  B.  EFFECTS  OF  CORRECTIONS 

Some  selected  results  of  applying  the  method  are  given  in  Figures  Bl 
through  B4;  two  solutions  for  each  of  the  four  days,  top  and  side  views  for 
each.  Original  track  from  the  first  member  of  the  array  pair  is  marked  with  a 
cross,  and  from  the  second  member  with  a  small  circle.  The  (timing) 
corrected  tracks  are  marked  with  connected  solid  lines  for  the  former  and 
dashed  lines  for  the  latter. 

The  first  set  in  Figure  Bl  shows  no  noticeable  corrections,  5  is  small.  The 
second  set  in  that  figure  has  8  =  3.54  and  is  a  real  data  version  of  the  situation 
depicted  by  Figure  3.  The  first  set  of  Figure  B2  appears  as  if  one  track  is 
corrected  more  than  the  other.  But  this  can  be  explained  because  the 
''stretching  angles"  are  different.  The  second  set  of  this  figure  shows  a  rather 
common  condition  in  that  the  corrected  track  is  at  a  shallower  depth  than  the 
original. 

The  first  set  in  Figure  B3  has  something  of  a  "showcase"  flavor;  the  top 
view  has  corrected  track  with  desirable  coherence.  (But  a  systematic 
separation  in  the  vertical  remains.)  The  second  set  also  shows  tantalizing 
improvement.  The  first  set  in  Figure  B4  has  8  =  -0.8,  and  provides  an 
example  of  modest  "contracting"  of  points.  The  second  set  is  an  extreme  case 
of  curved  track.  The  top  view  has  an  excellent  correction  displayed,  but  the 
side  shows  that  the  vertical  is  still  separated. 


30 


-400 

■5  -410 

H. 

S   -420 


-430 


5.26 


Top  View 


T  =  22 


Down  Range 

8  =  0.10 

Side  View 


Pair  (5,6) 


5.28 


5.3  5.32  5.34 

Down  Range 


5.36 


5.38 
xl04 


•1780 


S-1800- 

1-1820- 

U 

-1840- 


■380 


-400- 


Q   -420 


Top  View 


2.26  2.28  2.3  2.32  2.34 

Down  Range 
T  =  37  5  =  3.54  Pair  (1,2) 

Side  View 


2.36 


2.38 

X104 


**^\*- - 


-6         O    O  « 


>>' *  *>i  **f>*  >*  *  i    * 


«     w  \* 


x        x1*  * 


0\  o  e 


*c 


-440 


2.26  2.28 


2.3  2.32  2.34 

Down  Range 


2.36 


2.38 

Xl04 


Figure  Bl.   Before  and  After  Offset  Correction  -  3/23/89 


31 


-1880 
u 

£-1900 

g-1920 

u 


-1940 


2.24 


•360 


-370- 
&  -380 


■5 


-390 


2.24 


•  0^*. 


Top  View 


•  "  o 


•  a   O 


*        «.  1- 


XX    * 


X  XX  x 


X    XX    X  X   XX 


X  X   x   XX   x 


X  x  xxxxxx 


XXXXXxXxXx*    XXXXXi 


2.26 


T=43 


2.28 


2.3 

2.32              2.34 

2.36 

Down  Range 

xl04 

6=3.04 

Pair  (1,11) 

Side  View 

lt*A 


2.26 


2.28 


2.3 
Down  Range 


2.32 


2.34 


2.36 

Xl04 


3500 


2500 


-360 


-370 


r 


i3 


-340 


Top  View 


5.05      5.06      5.07      5.08      5.09       5.1       5.11 

Down  Range 


T=50 


6=4.91 
Side  View 


5.12      5.13      5.14      5.15 

Xl04 

Pair  (5,56) 


^^i^'t-^^s^^^^:^^^^- 


*\& 


5.05      5.06      5.07      5.08      5.09       5.1       5.11      5.12      5.13      5.14      5.15 

Down  Range  xlO 


Figure  B2.   Before  and  After  Offset  Correction  -  5/10/89 


-2800 
u 

c 

<*-2850 


e  c    ©  o  o 


Top  View 


e  o  e  e  o 


*X*X   XXX 


o  ©   e  © 


*  Xx* 


••<•«• 


u 


-2900 


2.46 


2.48 


T=50 


x-    X  X  X   X  x 


2.5 


•  *  •  ©    O 


•  ©  ©  © 


X  X  X  x  X   X 


•  •    6    ©  © 


X  X  X 


X    XX   X 


*  X  X  X 


x  x  x 


X  X  X 


2.52 

2.54              2.56 

2.58 

Down  Range 

X104 

6=2.72 

Pair  (2,11) 

-350 


Side  View 


-370 


2.46 


2.48 


2.5 


2.52 
Down  Range 


2.54 


2.56 


2.58 
xl(M 


0 
§   -500 

Qi 

1-1000 

u 

-1500 


-360 


£   -370 

Q. 

Q  -380 


-390 


Top  View 


«MHbM<**"*~*" 


X    X 


3.815      3.82      3.825      3.83 


T=50 


3.835  3.84 
Down  Range 
6=4.03 

Side  View 


3.845      3.85 


Pair  (3,4) 


3.855      3.86 

X104 


a*.    »    .  ° 

a^  ©         *  ft  "   °      ^ 

•  *  ~x!       °    x         -  * 


*    X* 


x      x      „ 


3.815       3.82       3.825      3.83      3.835       3.84      3.845      3.85       3.855       3.86 

Down  Range  xl(M 

Figure  B3.   Before  and  After  Offset  Correction  -  6/6/89 


33 


-390 


•5   -4W 

H. 


a . 


410 


-420 


Top  View 


3.86         3.88 


3.9  3.92         3.94         3.96         3.98 

Down  Range 

T=50    6=-0.80  Pair  (4,54) 

Side  View 


4.02 

X104 


WW    w         ^ 


*      X       X  ¥V  *     v>      *      v       >A*, 


3.86         3.88  3.9  3.92         3.94         3.96         3.98  4  4.02 

Down  Range  xlCM 


-5500 


1-6000 
'  -6500 


c 

u 


-7000 


5.705 


-340 


5   -350 

E. 


-360 


-370 


Top  View 


*> 


«• 


»•« 


* 


*S 


*  °«a 


** 


N^V 


x* 


?** 


**         X 


5.706 
T=50 


5.706  5.707  5.707 

Down  Range 
8=2.24 


5.708 
Pair  (15,16) 


5.709 
xl(M 


Side  View 


->"te\<«' 


* 


** V*  f***^ 


5.705  5.706  5.706  5.707  5.707 

Down  Range 


5.708 


5.709 

Xl04 


Figure  B4.   Before  and  After  Offset  Correction  -  7/21  /89 


34 


APPENDIX  C.  COVARIANCE  MATRICES 

The  covariances  of  the  residual  processes,  eq.  (38),  are  recorded  here  along 
with  M,  eq  (39),  which  estimates  M,  eq.  (16).  The  cross  covariances  Dxy  are 
converted  to  correlations,  Rxy,  by 

The  order  of  presentation  is  that  o(  Table  1.  The  number  of  points  used,  T,  is 
occasionally  smaller  than  its  Table  1  counterpart.  This  is  due  to  outlier 
rejection  based  upon  a  visual  view  of  the  track.  In  three  instances  these 
values  are  zero  and  then,  the  covariances  are  not  computed. 

In  a  number  of  instances  the  covariances  are  marked  with  a  "c"  or  "cc"  in 
the  identifying  column.  This  means  the  track  was  curved  (c)  or  excessively 
cu-ved  (cc)  so  that  the  straight  line  filter  could  not  reasonably  be  assumed  tc 
provide  a  valid  representation  of  the  path.  It  is  a  curiosity  that,  in  many  of 
these  instances,  the  cross  correlations  are  strong  and  their  use  produces,  via 
eq.  (16),  compensatory  values.  That  is,  the  M  matrices  do  not  appear  to  be 
especially  large.  Some  large  cross  correlations  also  appear  in  cases  that  pass 
the  visual  test  for  straight  line  track. 
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Table  C.l.   Covariances  and  Cross  Correlations  of  Residuals 
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APPENDIX  D.  AUTOCORRELATIONS 

Examples  oi  autocorrelation  functions  of  residuals  R(h)/R(0),  (see  eq.  (41)) 
are  presented  in  Figure  Dl  through  D6.  The  first  two  are  for  the  displacement 
process,  eq.  (40).   The  latter  four  are  for  the  three  components  of  displacement. 
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Figure  Dl.  Autocorrelation  of  noise  displacement 
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Figure  D3.    Autocorrelation  of  noise  components 
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Figure  D4.  Autocorrelation  of  noise  components 
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Figure  D5.    Autocorrelation  of  noise  components 
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APPENDIX  E.  COMPUTER  SOURCE  CODES 

The  major  programs  needed  to  produce  Table  1  are  documented  here. 

Basically  there  are  three  steps: 

(i)      Extract  pertinent  crossover  data  from  T-files 

(ii)    Produce  the  covariance  matrices  of  residuals;  DX/  Dy/  Dxy  and  M.  See 
equations  (38)  and  (39) 

(iii)  Develop   the  estimates,  eq.   (9)   and   (10).      Develop   supporting 
statistics,  eq  (17),  (18),  (19) 

Also  the  user  needs  access  to  Table  E.l,  the  coordinates  of  the  various 
arrays  at  the  Nanoose  range. 

Step  (i)  is  managed  by  the  program  KEYGATE,  written  in  FORTRAN  77, 
Miscosoft  optimizing  compiler  4.01.  This  program  reads  the  NUWES  T-files 
and  performs  a  series  of  gating  operations  in  order  to  produce  crossover  data 
in  the  required  seven  column  format.    It  will  prompt  the  user  for 

(a)  name  oi  the  range  (e.g.,  Nanoose) 

(b)  Number  of  records   to   ignore   (e.g.,   31   for  bypassing   the  DVT 
information) 

(c)  The  target  vehicle  mode  (e.g.  7) 

First,  it  will  strip  out  all  data  other  than  mode  2  or  mode  7.  Second,  it 
identifies  all  point  counts  for  which  a  position  vector  is  available  for  two  or 
more  arrays.  Third,  it  reads  an  array  location  file  and  removes  data  that 
cannot  be  reasonably  identified  with  an  overlap  region.  Next,  it  organizes  the 
data  by  pairs  oi  arrays.  Finally,  it  prompts  the  user  for  his  array  pair  and 
output  filename;  it  selects  and  records  all  the  data  of  the  desired  type. 
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Step  (ii)  is  contained  in  the  program  PRINCOM3.M.  This  is  a 
PCMATLAB  program  that  performs  the  eigen  analysis,  eq.  (32);  develops  the 
smoothed  fit,  eq.  (37);  and  the  covariance  matrices  (38)  and  (39).  The  latter  is 
placed  on  an  output  file  to  be  used  in  Step  (iii).  It  also  develops  the 
autocorrelation  sequence,  e.g.  (41)  of  the  displacement  process. 

A  slight  modification  of  this  program  can  be  made  to  develop 
autocorrelations  for  the  three  components  of  residuals  instead  of  for  the 
displacement  process. 

Step  (iii)  is  performed  by  the  program  TIMECOR,  a  FORTRAN  77  code. 
The  Microsoft  FORTRAN  optimizing  compiler  4.01  is  used.  The  user  should 
note  that  requirements  of  the  three  input  files: 

VELOCITY.DAT  is  a  two-column  data  set  containing  the  water  layer 
boundaries  and  the  sound  speeds  in  25  ft  increments 

TRPDOTRX.DAT  is  the  seven-column  crossover  data  output  of  KEYGATE 
to  which  has  been  appended  the  locations  of  the  two  arravs  as  the  last 
record. 

MMATRIX.M  is  the  covariance  matrix,  M,  produced  by  PRINCOM.M  in 
Step  (ii). 

Since  the  ray  tracing  exit  layer  elevation  angles  are  not  contained  in 
T-files  it  is  necessary  to  reconstruct  them  in  order  to  compute  Equation  (5). 
This  is  accomplished  by  the  subroutine  TGEN.  The  methodology  of  TGEN  is 
explained  in  [3]  under  the  heading  of  ray  fitting. 

There  are  two  output  files  which,  taken  together,  contain  all  of  the 
information  indicated  in  Table  1. 
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TABLE  E.l.  COORDINATES  OF  THE  NANOOSE  ARRAYS 


Array  Number 

Date  of  Survey 

Downrange 

Crossrange 

Depth 

0 

6/20/85 

12188.01 

-131.52 

-1295.33 

1 

6/20/85 

19463.16 

-174.99 

-1308.76 

2 

7/12/85 

26991.39 

-109.83 

-1323.25 

3 

1/07/88 

34505.10 

-80.76 

-1323.32 

4 

10/24/88 

42005.19 

-55.17 

-1318.28 

5 

6/20/85 

49497.00 

-25.23 

-1315.58 

6 

6/20/85 

56972.28 

-21.21 

-1308.50 

7 

7/30/85 

64680.66 

15.33 

-1353.39 

8 

11/16/88 

71969.73 

-29.28 

-1300.89 

9 

5/07/84 

3.00 

3.00 

1.00 

10 

3/12/84 

47100.00 

-3600.00 

-1300.00 

11 

7/18/85 

23173.89 

-6488.40 

-1312.09 

12 

6/20/85 

30731.25 

-6553.05 

-1312.90 

13 

6/20/85 

38213.61 

-6640.77 

-1323.05 

14 

6/20/85 

45647.07 

-6513.18 

-1324.78 

15 

6/19/85 

53249.43 

-6354.60 

-1316.66 

16 

9/13/85 

60859.74 

-6356.07 

-1313.42 

17 

6/16/87 

68217.93 

-6524.10 

-1313.43 

54 

2/02/88 

38029.95 

5401.98 

-1212.69 

55 

6/20/85 

45645.75 

6369.66 

-1188.12 

56 

7/30/85 

53180.13 

6417.96 

-1218.84 

57 

7/30/85 

60745.71 

6419.40 

-1088.24 

23 

6/20/85 

41605.14 

-12150.18 

-1268.23 

24 

4/17/89 

49572.00 

-12966.00 

-1300.00 

25 

10/24/88 

56993.79 

-12999.33 

-1205.48 

26 

8/08/88 

64442.94 

-12971.04 

-1255.35 

27 

7/15/80 

22119.60 

-15908.70 

83.00 

28 

5/04/83 

45000.00 

1500.00 

-1350.00 

29 

2/02/79 

.00 

.00 

.00 

KEYGATE 

PROGRAM  KEYGATE 


Program  to  read  in  raw  data  from  Keyport  hydrophone  arrays,  segregate  it  by  mode,  and  throw  out 
unusable  records.  The  output  of  this  program  is  to  be  read  in  by  the  program  KEYMAIN. 

Modified  by  Colin  R.  Cooper      11/15/89 

1NTEGERM  CRT,  KBD 
CHARACTER*25  DSNAME,  SITNAM 
CHARACTERS  TEMPI,  TEMP2,  TEMP3,  TEMP4 

PARAMETER  (KBD=5,  CRT =6) 

WRITE(CRT,»)  '  Please  enter  the  name  of  your  input  file:  ' 
REAEKKBD/(A)')  DSNAME 

WRITE(CRT,*)  '  Please  enter  the  name  of  the  range  configuration  file:', 
READ(KBD,'(A)')  SITNAM 

TEMPI  =  TEMPI  TMP 
TEMP2  =  TEMP2TMP 
TEMP3  =  TEMP3.TMP 
TEMP4  =  TEMP4.TMP 

CALL  STRMOD(CRT,KBD,TEMPl,DSNAME) 

CALL  PAIR(CRT,KBD,TEMP1,TEMP2) 

CALL  RANGE(CRT,KBD,SITNAM,TEMP2,TEMP3) 

CALL  PAIR2(CRT,KBD,TEMP3,TEMP4) 

CALL  REC(CRT,KBD,TEMP4) 

VVRITEtCRT,*)'  Operation  complete     KEYGATE  terminating...  ' 

STOP 
END 


SUBROUTINE  STRMOD(CRT,KBD,TEMPl, DSNAME) 


Program  to  strip  all  modes  except  2's  and  7's  from  the  Keyport  data.  (2  indicates  target  ship,  7 
indicates  torpedo). 


CHARACTER  DO*2,  DSNAME*25,  TEMPI *9 
INTEGER  PC,  ARRAY,  CRT,  NHEAD 

OPEN(l,FlLE=DSNAME,STATL'S=,OLD) 

WRITE(CRT,*)'  How  many  records  of  header  do  you  want  to  strip  off  the  file7', 
READ(KBD,*)  NHEAD 
DO  11  1  =  1,  NHEAD 
READd,*) 
11     CONTINUE 
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WRITEfCRT,*)'  Input  mode  to  be  kept?' 
READ(KBD,*)  NUM 

OPEN(2,FILE=TEMPl/STATUS=,NEW) 

10     READn,100,END=50,ERR=40)PC/DO,X,Y,Z,ARRAY,MODE 

IRDO.NE. '  ,)GOTO20 

IF(MODE  .NE.  NUM)  GOTO  20 

YVRITE(2,1 1 0)PC,X,Y,Z,  ARRAY,MODE 
20    CONTINUE 

GOTO  10 
40     WRITE(CRT,*)'  THERE  WAS  AN  ERROR  IN  THE  FILE' 
50    CONTINUE 

100     FORMATd5,A2,lX,F7.1,2X,F7.1,2X,F7.1,30X,I2,2X,Il) 
110     FORM  AT(1X,I5,2X,3F1 0.1, 2X,I2,2X,I2) 

CLOSE  (UNIT=1) 

CLOSEdJNIT=2) 

RETURN 

END 


SUBROUTINE  PAIR(CRT,KBD,TEMP1,TEMP2) 


Program  to  pair  point  counst  after  the  data  has  been  gated  by  STRMOD.  Second  pass. 


DIMENSION  X(200),  Y(200),  Z(200) 

INTEGERM  PC(200),  ARRAY(200),  MODE(200),  CRT,  HOLD 

CHARACTERS  TEMPI,  TEMP2 

OPENd  ,FILE=TEMPl,STATUS='OLD') 
OPEN(2,FILE=TEMP2,STATUS=,NEW) 

HOLD  =  0 
IREC  =  0 
NREC  =  0 
1  =  1 


Read  records  by  two's  to  compare  point  counts. 


READ(1,\END=40,ERR=30)  PC(I),  X(l),  Y(I),  Z(I), 
+  ARRAY(I),  MODE(I) 

10     READd,',END=40,ERR=30)  PCd+1),  X(I+1),  Y(I+1),  Z(I+1), 
+  ARRAYd+1),  MODEd+1) 

NREC  =  NREC  +  1 

IF(PC(I+1)  .EQ.  HOLD)  THEN 

WRITE(2,100)  PCG+1),  Xd+1),  Y(I+1),  Z(I+1), 
+  ARRAYd+1),  MODEd+1) 

HOLD  =  PC(l+l) 

GO  TO  20 
END  IF 

IF(PC(I).EQ.  PC(I+1))THEN 

WRITE(2,100)  PC(I),  X(I).  Y(I),  Z(I),  ARRAY(I),  MODE(I) 

v\T<rTE(2,ioo)  pcd+u,  xd-i),  Yd+i),  za+i), 
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+  ARRAYU+1),  MODEO+1) 

HOLD  =  PC(I+l) 
IREC  =  IREC  +  1 
GO  TO  20 
END  IF 

IF(PC(I)  .NE.  PC(I+D)  THEN 

pcd)  =  pca+i) 

xd)  =  xa+i) 
Y(i)=Ya+i) 

za)  =  za+i) 

ARRAY(l)  =  ARRAY0+1) 
MODE(l)  *  MODE(I+l) 
1  =  1 
END  IF 

GOTO  10 
20    1  =  1  +  1 
GO  TO  10 

30     WRITE(CRT,*)'  THERE  IS  AN  ERROR  IN  THE  DATA  RLE  IN  REC.',NRE 
WRITE(CRT,»)'  ...  OPERATION  TERMINATING  DUE  TO  ERROR.  ' 
STOP 

40    CONTINUE 

CLOSE(UNIT=l) 
CLOSE(UNIT=2) 

100    FORMAT(lX,I5,2X,F7.1,2X,2F9.1,2X,I2,2X,I2) 

RETURN 

END 


SUBROUTINE  RANCE(CRT,KBD,SITNAM,TEMP2,TEMP3) 

This  program  completes  the  third  gating  of  Keyport  range  data.  It  reads  array  location  data  from  a  site 
specific  configuration  file  and  tests  to  see  if  the  data  is  in  the  valid  overlap  area. 


INTEGERS  ARRAY,  CARRAY,  CRT,  P 
REALM  CONFIG(200,4),  LX,  LY,  LZ,  MAXVAL 
CHARACTER  SITNAM*25,  TEMP2*9,  TEMP3*9 


Oper.  input  and  output  files: 


OPENd  ,FlLE=TEMP2,STATUS='OLD) 

OPENC2,FILE=SITNAM/STATUS=OLD') 

OPEN(3,nLE=TEMP3,STATUS=,NEV\") 


Read  site  configuration  into  CONFIG  array: 


NREC  =  0 
1  =  1 
10     READ(2,*,END=40,ERR=30)  CONFlG(l,1 ),  CONFlC(I,2),  CONFIG(I,3), 
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+  CONFlC(I,4) 

NREC  =  NREC  +  1 

1  =  1  +  1 

GOTO  10 
30     W'RITE(CRT,*)'  There  was  an  error  reading  the  config  file  in  record',NRB 
40    CONTINUE 

CLOSE(UNIT=2) 

NDREC  =  0 

Read  X,  Y,  Z,  and  ARRAY  from  input  data  file: 

45     READa,\ERR=70,END=80)  PC,  X,  Y,  Z,  ARRAY,  MODE 
DO50I  =  l,NRE 
LX  =  CONFlG(I,l) 
LY  =  CONFIG(I,2) 
LZ  =  CONF1G(I,3) 
CARRAY  =  INT(CONFlG(I,4)) 
MAXVAL  =  4700. 

Match  array  number  in  data  file  with  that  in  config.  file.   If  they  are  equal  compute  slant  range 
distance  (SRDIST): 

IRARRAY  .EQ.  CARRAY)  THEN 

SRDIST  =  SQRT(  (X  -  LX)~2  +  (Y  -  LY)"2  +  (Z  -  LZ)*»2  ) 
IF(MAXVAL  .GE.  SRDIST)  THEN 

WRJTE(3,100)  PC,  X,  Y,  Z,  ARRAY,  MODE,  SRDIST 
NDREC  =  NDREC  +  1 
END  IF 
END  IF 

50    CONTINUE 
GOTO  45 

70     WRrTE(CRT,*)'  There  is  an  error  in  the  data  file.  ' 
80    CONTINUE 

CLOSE(UNTT=l) 
CLOSE(UNTT=3) 


100    FORMAT(3X,I5,3(3X,F10.U2(3X,I2UX,F8.2) 

RETURN 
END 


SUBROUTINE  PAIR2(CRT,KBD,TEMP3,TEMP4) 


Program  to  pair  point  counts  after  the  data  has  been  tested  by  RANGE.   Fourth  pass. 


DIMENSION  X(200),  Y(200),  Z(200),  SRDIST(200) 
INTEGERM  PC(200),  ARRAY(200),  MODE(200),  CRT,  HOLD 
CHARACTERS  TEMP3,  TEMP4 
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OPENO  ,FILE=TEMP3,STATUS='OLD') 
OPEN(2/nLE=TEMP4,STATUS=,NEW) 

HOLD  =  0 
IREC  =  0 
NREC  =  0 
1  =  1 


Read  records  by  two's  to  compare  point  counts. 


READ(1,*,END=40,ERR=30)  PC(I),  X(I),  Y(I),  Z(I), 
+  ARRAYd),  MODEd),  SRDIST(I) 

10     READ<1,*,END=40,ERR=30)  PCd+1),  X(I+1),  Yd+1),  Z(I+1), 
+  ARRAY(I+1),  MODEd+1),  SRDIST(l-t-l) 

NREC  =  NREC  +  1 

IF(PC(I+1)  .EQ.  HOLD)  THEN 

WRITE(2,100)  PCd+1),  X(I+1),  Y(I+1),  Z(I+1), 
+  ARRAYd  +  1),  MODEd+1),  SRDlSTd+1) 

HOLD  =  PCd+1) 

GO  TO  20 
END  IF 

IF(PC(I).EQ.PC(I+1))THEN 

WRITE(2,100)  PCd),  X(I),  Y(I),  Z(I),  ARRAYd),  MODEd), 
+  RDIST(I) 

WRITE(2,100)  PCd+1),  Xd+1),  Y(I+1),  Zd+1), 
+  ARRAYd+1),  MODEd+1),  SRDIST(I+1) 

HOLD  =  PCd+1) 

IREC  =  IREC  +  1 

GO  TO  20 
END  IF 

IF(PC(I)  .NE.  PC(1+1))THEN 

PC(1)  =  PC(1+1) 

X(l)  =  Xd+l) 

Y(1)  =  Y(I+1) 

Z(1)  =  Z(I+1) 

ARRAY(l)  =  ARRAYd+1) 

MODE(l)  =  MODEd+1) 

SRDIST(l)  =  SRDIST(I  +  1) 

1  =  1 
END  IF 

GOTO  10 
20     1  =  1+1 
GO  TO  10 

30     WRrTE(CRT,»)'  THERE  IS  AN  ERROR  IN  THE  DATA  FILE  IN  REC.,NRE 
WRITE(CRT,')   ...  OPERATION  TERMINATING  DUE  TO  ERROR.  ' 
STOP 

40     CONTINUE 

CLOSE(UNTT=l) 
CLOSE(UNTT=2) 

100    FORMAT(lX,I5,2X/F7.1,2X,2F9.1,2X/I2/2X,12,2X/F8.2) 
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RETURN 
END 


SUBROUTINE  REC(CRT,KBD,TEMP4) 


Program  to  produce  the  final  gating  of  Keyport   hydrophone  array  test  data. 


INTEGER  PC,  ARRAY,  PCI,  PC2,  Al,  A2,  ARRAY1,  ARRAY2,  CRT 
INTEGER  PCA(IO),  ARRAYA(IO),  HOLD 
DIMENSION  XAOO),  YAOO),  ZA(lO),  SRDISTAOO) 
CHARACTER     OUTFIL'25,  ANSM,  TEMP4*9,  TEMPI *9,RANCE»7 

RANGE='NANOOSE' 

99     WRITE(CRT,*)'  What  is  the  name  you  wish  to  give  to  ', 
+  '  your  output  file? 

READ(CRT,'(A)')  OUTFIL 

WRITE(CRT,*)'  Enter  the  numbers  of  the  arrays  to  be  paired:  ' 
READ(KBD')  Al,  A2 

OPENd  ^LE^EMP^STATUS^OLD1) 
OPEN(2,FILE=TEMPl.DAr/STATUS=,OLD') 

5    READ(1,»,END=8,ERR=8)  PC,  X,  Y,  Z,  ARRAY,  MODE,  SRDIST 
IF(( ARRAY  .EQ.  Al).OR.( ARRAY  .EQ.  A2))  THEN 

WRITE(2,100)  PC,  X,  Y,  Z,  ARRAY,SRDIST 
END  IF 
GOTO  5 

8    CONTINUE 
CLOSE  (UNIT=1) 


Routine  to  pair  data  again: 


REWIND  2 
OPEN(4,FILE=TEMP2.DAT',STATUS='OLD) 

1  =  1 

IFLAG  =  0 
FIRST  =  1 
HOLD  =  0 

M  =  0 

9    READ(2,100,END=40)  PCA(I),  XA(I),  YA(I),  ZA(I), 
+  ARRAYA(I),  SRDISTA(I) 

IF(RRST.EQ.  DTHEN 
FIRST  =  0 
HOLD  =  PCA(I) 
END  IF 

IF(PCA(I)  .EQ.  HOLD)  THEN 
1  =  1  +  1 

M  =  M  +  1 
GOTO  9 


ELSE 


In  cases  where  there  are  three  or  more  reports  for  a  given  point  count,  segregate  by  comparing 
SRDIST.   if  there  are  more  than  3  reports,  discard  all. 

IF  (M  .EQ.  3)  THEN 

IF(ARRAYA(M)  .EQ.  ARRAY A(M-l))  THEN 
IF  (ARRAYA(M)  .EQ.  ARRAYA(M-2))  THEN 
GOTO  50 
ELSE 

IF  (ABS(SRDISTA(M-2)-SRDISTA(M))  .LT. 
+  ABS(SRDISTA(M-2)-SRDISTA(M-l)))  THEN 

WRrTE(4,100)  PCA(M),  XA(M),  YA(M),  ZA(M), 
+  ARRAYA(M),  SRDISTA(M) 

WRITE(4,100)  PCA(M-2),  XA(M-2),  YA(M-2),  ZA(M-2), 
+  ARRAYACM-2),  SRDISTA(M-2) 

IFLAG  =  1 
ELSE 

WRITE(4,100)  PCA(M-l),  XA(M-l),  YA(M-l),  ZA(M-l), 
+  ARRAYA(M-l),  SRDISTA(M-l) 

\VRITE(4,100)  PCA(M-2),  XA(M-2),  YA(M-2),  ZA(M-2), 
+  ARRAYACM-2),  SRDISTA(M-2) 

IFLAG  =  1 
END  IF 
END  IF 
ELSE 

IF  (ARRAYA(M)  .EQ.  ARRAYACM-2))  THEN 
IF  (ABS(SRDISTA(M-l)-SRDISTA(M))  .LT. 
+  ABS(SRDISTA(M-l)-SRDISTA(M-2)))  THEN 

WRITER,  100)  PCA(M),  XA(M),  YA(M),  ZA(M), 
+  ARRAYA(M),  SRDISTA(M) 

WRITE(4,100)  PCA(M-l),  XA(M-l),  YA(M-l),  ZA(M-l), 
+  ARRAYA(M-l),  SRDISTA(M-l) 

IFLAG  =  1 
ELSE 

WRITE(4,:iOO)  PCA(M-l),  XA(M-l),  YA(M-l),  ZA(M-l), 
+  ARRAYA(M-l),  SRDISTA(M-l) 

VNrRITE(4,100)  PCACM-2),  XA(M-2),  YA(M-2),  ZA(M-2), 
+  ARRAYA(M-2),  SRDlSTA(M-2) 

IFLAG  =  1 
END  IF 
ELSE 

IF  (ABS(SRDISTA(M)-SRDISTA(M-D)  .LT. 

ABS(SRDISTA(M)-SRDISTA(M-2)))  THEN 
VVRITEi4,lO0)  PCA(M),  XA(M),  YA(M),  ZA(M), 
+  ARRAYA(M),  SRDISTA(M) 

WRITE(4,100)  PCA(M-l),  XA(M-l),  YA(M-l),  ZA(M-l), 
+  ARRAYA(M-l),  SRDISTA(M-l) 

IFLAG  =  1 
ELSE 

\VRITE(4,100)  PCA(M),  XA(M;,  YA(M),  ZA(M), 

ARRAYA(M),  SRDISTA(M) 
WRITER,!  00)  PCA(M-2),  XA(M-2),  YACM-2),  ZA(M-2), 
+  ARRAYA(M-2),  SRDISTA(M-2) 

IFLAC  =  1 
END  IF 
END  IF 


END  IF 


ELSE 

IF  ((M  .EQ.  2,  .AND.  (ARRAYA(M)  .NE.  ARRAYA(M-l)))  THEN 
WRITE(4,100)  PCA(M),  XA(M),  YA(M),  ZA(M),ARRAYA(M), 
+  SRDISTA(M) 

WRITE(4,100)  PCA(M-l),  XA(M-l),  YA(M-l),  ZA(M-l), 
+  ARRAYA(M-l),  SRDISTA(M-l) 

END  IF 
END  IF 
50  CONTINUE 

25  CONTINUE 

IFLAG  =  0 

PCA(1)  =  PCA(I) 

XA(1)  =  XA(I) 

YA(1)  =  YA(I) 

ZA(1)  =  ZA(I) 

ARRAYA(l)  =  ARRAYA(I) 

SRDISTA(l)  =  SRDISTA(I) 

1  =  2 

M  =  l 

HOLD=PCA(l) 

DO  45  L  =  2,10 
PCA(L)  =  0 
XA(L)  =  0 
YA(L)  =  0 
ZA(L)  =  0 
ARRAYA(L)  =  0 
SRDISTA(L)  =  0 
45         CONTINUE 

END  IF 

GOTO  9 
40    CONTINUE 

NREC  =  0 

CLOSE  (UNIT=4) 

OPEN(7/nLE=,TEMP2.DAT/STATUS=,OLD) 
OPEN(9,FILE=OUTFlL,STATUS='NEW) 

WRITE(9,300)  RANGE,  Al,  A2 

Read  in  array  data  in  two  record  pairs: 

10    READ(7,100,END=60,ERR=70)  PCI,  XI,  Yl,  Zl,  ARRAY1,  SRDISTl 
READ(7,100,END=60,ERR=70)  PC2,  X2,  Y2,  Z2,  ARRAY2,  SRDIST2 

IF(ARRAY1  .EQ.  Al  .AND.  ARRAY2  .EQ.  A2)  THEN 


If  arrays  are  in  specified  order  (e.g.  4,5): 


WRITE(9,200)  PCI,  XI.  Yl,  Zl,  X2,  Y2,  Z2 
END  IF 

IF(ARRAY1  .EQ.  A2  .AND.  ARRAY2  .EQ.  AT)  THEN 


If  array's  are  in  reverse  order  (eg.  5,4): 


WRrrE(9,200)  PCI,  X2,  Y2,  72,  XI,  Yl,  Zl 
END  IF 


Increment  record  counter: 


NREC  =  NREC  +  1 

GOTO  10 

70     VVRITEtCRT,*)'  There  is  a  bad  record  in  the  file.' 
60    CONTINUE 

CLOSE(UNIT=2) 
CLOSE(UNIT=7) 
CLOSE(UNTT=9) 

100    FORMAT(I5,3(2X,F8.1),2X,I2,2X,F8.2) 
200    FORMAT(2X,I5,lX,6(2X,Fll.D) 
300     FOR\.AT(16X,A10,2X,I2,3X,I2) 

VVRITEtCRT*)  '  Do  you  want  to  try  another  array  pair?  (Y/N)' 

READ(KBD,'(A)')  ANS 

IRAN'S  .EQ  "Y  .OR.   ANS  .EQ.  V)  GO  TO  99 

RETURN 

END 


PRINCOM3 

function!  Au  toco  rx,Au  toco  rvl  =  princom3(f  name) 


3/13/90 
PRINCOM3.M    will  evaluate  the  principil  components  of  the  passed  data  file.   Computes  the 
autocorrelation  of  the  displacement  process. 


iname  =  ['d:\mdataV  fname  '.out']; 

evalU'load  '  iname  ]), 
evaKI'keyout  =  '  fname  ';']); 
eval([  clear  '  fnamel); 


^ 


X=keyout(:,5:7); 
Y=keyout(:,2:4); 

[len,n]=sLze(keyout); 


Compute  averages  and  covariance  matices  for  the  tracks  produced  by  each  of  the  arrays. 

Xav=ones(len,1)*sum(X)/len  ; 

Yav=onesflen,l)*sum(Y)/len  ; 

XX=X-Xav ; 
YY=Y-Yav ; 


ipc=keyout(:,l); 

CX=cov(XX); 
CY=cov(YY); 


Develop  the  eigen  analysis.   Choose  largest  eigen  value. 


|WX,DX)=eig(CX)  ; 

[WY,DY]=eig(CY)  ; 

ifDXO,l)>DX(2,2) 

ifDX(l,l)>DX(3,3) 
lx=l  ; 

end 
elseifDX(2,2)>DX(3,3) 

lx=2      ; 
end 
ifDX(3,3)>DX(l,l) 

if  DX(3,3)>DX(2,2) 
lx=3; 

end 
end 
ifDY(U)>DY(2,2) 

ifDY(l/l)>DY(3,3) 
ly=i ; 

end 
elseif  DY(2,2)>DY(3,3) 

ly=2      ; 
end 
ifDY(3,3)>DY(l,l) 

if  DY(3,3)>DY(2,2) 
lv=3; 

end 
end 

PX=XX*WX; 
PY=YY*WY; 

ux=PX(:,lx); 
uy=PY(:,ly); 


Modify  data  projection  onto  the  first  principal  component  to  account  for  constant  speed. 


ipcav=sum(ipc)/length(ipc); 


mx=(sum(ux.*(ipc-ipcav)))/(sum((ipc-ipcav).*(ipc-ipcav)))  ; 
rny=(sum(uy.*(ipc-ipcav)))/(sum((ipc-ipcav).*(ipc-ipcav)))  ; 

ax=(sum(ux)/length(ux))  -  (mx'ipcav), 
ay=(sum(uy)/length(uy))  -  (my*ipcav); 

uhx=(mx*ipc)+ax; 
uhy=(my*ipc)+ay; 

PPX=zeros(len,3); 
PPY=zeros(len,3); 

PPX(:,lx)=uhx; 
PPY(:,ly)=uhy; 


Compute  the  straight  line  tracks  in  the  original  coordinate  system. 


XSL=(PPX*WX')+Xav; 
YSL=(PPY*WY')+Yav; 


Develop  the  displacement  process  of  residuals. 


distx  =  sqrt((X(:,l)  -  XSL(:,1)).A2  +  (X(:,2)  -  XSL(:,2)).A2  +  (X(:,3)  -  XSL(:,3)).A2); 
disty  =  sqrt((Y(:,l)  -  YSL(:,1)).A2  +  (Y(:,2)  -  YSL(:,2)).A2  +  (Y(:,3)  -  YSL(:,3)).A2); 

davx  =  sum(distx)/length(distx); 
davy  =  sum(disty)/length(disty), 

\T(  =  distx  -  davTc; 
vy  =  disty  -  davy, 


Compute  the  auto  correlations. 


Autocorx  =  zeros(length(vx),l); 
Autocory  =  zeros(length(vy),l); 

for  k  =  l:length(vy) 

for  1  =  1  :length(vy)  -  k  +  1 

Autocorx(k)  =  Autocorx(k)  +  vx(i)*vx(i  +  k-1); 
Autocory(k)  =  AutocoryOO  +  vy(i)*vy(i  +  k-1); 
end 

Autocorx(k)  =  Autocorx(k)/sum(vx.A2); 
Autocory(k)  =  Autocory(k)/sum(vy.A2); 
end 

clg 

subplot(211) 

plot(Y(:,l),Y(:,2);or',YSL(:,l)/YSL(:,2)/-r'/X(;,l),X(:/2);xg',XSL(:,l),XSL(:/2);-g') 

title([fname  '  -  Corrected  Tracks  with  Principal  Components']) 

xlabcK'Down  Range'), ylabeKCross  Range*) 

subplot(212) 

plot{Y(:,l),Y(:,3);or,/YSL(:/l)/YSL(:/3)/,-r,,X(:,l),X(:/3);xg'J  XSL(:,l),XSL(:,3)/-g') 

titifr'Corrected  Tracks  with  Principal  Components') 

xlabcK'Down  Range'), yiabeK 'Depth') 
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pause 
clg 

h  =  0:length(Autocorx)-l; 

subp)ot(21 1  ),plot(h,Autocory,'or',h,Autocory,'-r') 

titled fname  '  -  Autocorrelation  of  Distance,  Y  Array' l),grid 

subplot(212),plot(h,Autocorx,'or',h,Autocorx,'-r') 

titlefAutocorrelation  of  Distance,  X  Array'),grid 

clg 

plot(h(l  :27),Autocorx(l  :27),'or',h(l  :27),Autocorx(l  :27),'-r') 

titleCAutocorrelation  of  Residual  Distance'), grid 

xlabelCPoint  Count  Lag') 


TIMCOR 

PROGRAM  TIMECOR 


06/06/90 
This  program  estimates  the  timing  synchronization  offset  and  drift  parameters  based  upon  cross-over 
data  from  T-files.  Since  exit  angles  and  transit  times  are  not  recorded  to  these  files  they  must  be  re- 
estimated.    Location  information  for  each  of  the  involved  arrays  is  also  required. 


This  program  was  compiled  using: 

Microsoft  FORTRAN  OPTIMIZING  COMPILER  ver.  4.01 


This  files  must  be  able  to  access  the  following  files 

VELOCITY.DAT      -  Sound  velocity  vs.  depth  data. 

TRPDOTRX.DAT    -  Torpedo  tracking  data. 

MMATRIX.M  -  Correlation  matrices  for  data. 

OUTPUTl.COV       -  Output  file  containing  delta  and  tnot 

OUTPUT2.COV       -  Output  file  containing  cross  correlation  and  M  matrices. 

Users  Notes: 

-  VELOCITY.DAT  file  contains  layer  boundaries  and  sound  velocities  for  25  ft.  depth  increments. 

-  TRPDOTRX.DAT  file  contains  point  counts  and  three  position  components  for  each  of  the  two 

contributing  arrays.   The  position  data  from  the  lower  numbered  array  is  columns  two  through 
four.  The  last  record  contains  the  coordinate  position  of  the  two  arrays.. 

-  MMATRIX.M  file  contains  the  set  of  covariance  matrices  produced  by  AUTOCOR3.M. 

-  The  two  output  files  contain  all  the  information  appearing  in  Table  1. 

DIMENSION  IPCO50) 

CHARACTERS  LINE 

REA L»8  Yl (1 50,2), Y2(l 50,2), Y3(l 50,2), YC1  (1 50,2), YC2(1 50,2) 

REAL*8  LL(55),G(55),VV(55),A2,Pl,P2,V0,V],TTIME(150,2) 

REAL*8THETA(150,2),DEPTH(55),PHI(150,2),V(150,2),DZ 

RE  AL*8Bl(150,2),B2n50,2),B3(l  50,2),  YC3(1 50,2) 

REAL*8  DEN,CG,SSB,SSR,TD,TM,SIC,DD,TEMP 

REAL*8  MEQ,GEQ(150),DELTA,T1,TNOT,D1,D2,MSR,MSB 

REAL*8  DUM1,DUM2,DUM3,DUM4,DUM5,DUM6,DUM7 

REAL*8  M(3,3),C(150,3),VDEL,SDEL,VM,SDVM,COVDM 
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REAL^DELTAN,MEQN,R,CX(3,3),CY(3,3),CXY(3,3),RXY(3,3)T 
INTEGERS  TEE 


Read  in  the  data  from  the  data  files. 


OPENCUNITs^nLE^VELOCITY.DAT^STATUS^OLD') 

OPEN(UNIT=7/FlLE=TRPDOTRX.DAT,/STATUS=,OLD) 

OPEN(UNnT=10/FILE='MMATRIX.M',STATUS=,OLD) 

OPEN(UNIT=n/nLE=,OUTPUTl.COV/STATUS=,OLD) 

OPEN(UNIT=12,nLE=,OUTPUT2.COV'/STATUS=,OLD') 

READ(10#120)CX(l#l)/CX(2,l),CX(3#l)fCX(l,2)/CX(2/2), 

CX(3,2),CX(1,3),CX(2,3),CX(3,3) 
READ(10,120)CY(l,l),CY(2,l),CY(3,l),O'(l,2),CY(2,2), 

CY(3,2),CY(1 ,3),CY(2,3),CY(3,3) 
REAEK10/120)CXY(l/l),CXY(2,l),CXY(3,l),CXAj'(l/2)/CXY(2,2)/ 

CXY(3,2),CXY(1,3)/CXY(2,3),CXY(3,3) 
REAEK10,124)T 
TEE  =  IDINTfT) 

IF  (TEE.EQ.O)  THEN 
0021  =  1,3 

RXY(1,I)  =  O.ODO 
RXY(2,I)  =  O.ODO 
RXY(3,I)  =  O.ODO 

2  CONTINUE 
ELSE 

D04I  =  1,3 
D03J  =  1,3 

RXY(J.D=CXY(J,I)/DSQRT(CX(I,I)»CY(JJ)) 
M(J.D  =  CXO.D  +  CY(J,I)  -  CXY(LJ)  -  CXY(J,D 

3  CONTINUE 

4  CONTINUE 
ENDIF 

5  READ(12,'(A)',END=7)LINE 
GOTO  5 

7  BACKSPACE  12 
WRITE(12,122)TEE 
DOS  I  =  1,3 

8  WRITE(12#121>CX(I,1),CX(I,2),CX(I,3),CY(I,1),CY(I,2). 

CY(l/3)/RXY(I,l),RXY(I/2),RXY(I,3),M(I,l)/M(I,2)/M(I,3) 
WRITE(12,123) 

1=1 

10     REALX7,')IPC(I;,Y1(I,1)/Y2(I,1),Y3(L1), 

IF(IPC(I).EQ.999)  GOTO  15 

1=1+1 

GOTO  10 
15     LEN=I-1 

120     FORMATOX,9(E15.4,1X)) 
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121  FORMAT(5X,I2,5X/l  ',25X/I  ',2Xt'\  ',25X,1  \2X;\  ',25X/I ' 

2X/I  \25X/I ') 

122  FORMAT(12X/'r',lX/3(F7.2,lX);T,2X/f,/lX,3(F7.2/lX),,l,/ 

2X/^lX,3(F7.2/lX)/r,2X,f,,lX,3(F7.2,lX)/T) 

123  FORMATCnX/L'^SX/j^X/L^X/J'^X/L'^SX/J' 

2Xri'f25X,'j) 

124  FORMAT(lX,El5.4) 


Read  the  VELOCITY.DAT  file  and  prepare  for  isogradient  raytracing. 


DZ=25 
1  =  1 
25     READ(2,',END=30)  TEMP,VV(I) 
1  =  1+1 
GOTO  25 
30    CONTINUE 

CC=(W(I-l)-VV(I-7))/6. 
DO  35  J  =  1,55 

VV(J)=VV(J-D+CG 
35    CONTINUE 
LL(1)=125 
DEPTH(1)=0 
DO  40  1=255 

LL(I)=LL(I-1)+DZ 
DEPTH(I)=DEPTH(I-1  )+DZ 
40    CONTINUE 
DO  43 1=2,55 

G(I-1)=(VV(I)-W(I-1))/DZ 
43    CONTINUE 

C(55)=C(54)+GG/DZ 
DO  45  1=152 
45    CONTINUE 
DO80J=U 
N=15 
IFLAC=0 
DO  70  IT=1,LEN 


Set  variables  the  the  call  to  TGEN  subroutine.   TGEN  returns  the  time  and  elevation. 

A2=Y3(LEN+1J) 

P1=SQRT((Y1(ITJ)-Y1(LEN+1/I))»*2+(Y2(IT,J)-Y2(LEN+1/J))*»2) 

P2=-Y3(IT,J) 

IF(IT.NE.l)  THETA(IT,J)=THETA(IT-1J) 

CALLTGEN(LL,G,VV,A2,Pl/P2/THETA(rr/J),TTIME(rT/J), 

V0,V1,DEPTH;IFLAG) 
IFLAG=1 


Find  the  azimuth  angles  from  each  arrav. 


PHI(IT,J)=DASIN((Y2(1TJ)-Y2(LEN+1J))/SQRT((Y2(1TJ)- 
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Y2(LEN+lJ))-2  +  (Y1(IT,J)-Y1(LEN+1,J))**2)) 
IF((Yl(rTJ)-Yl(LEN+lJ)).LT.O)  PHI (TTJ)=3.1 41 59265359 
-PHKITJ) 


Locate  the  layer  containing  the  source  and  set  it's  velocity. 


50  IF(-Y3(IT,J).LE.DEPTH(N).AND.-Y3(rTJ).GT.DEPTH(N-l))  GOTO  60 

IF(-Y30TJ).GT.DEPTH(N))  THEN 

N=N+1 
ELSE 
N=N-1 

IF(N.LE.l)  THEN 
N=l 

GOTO  60 
ENDIF 
END1F 
GOTO  50 
60  V(IT,J)=VV(N) 


Calculate  the  B(T)  adjustments   (Sperical  Coordinates  :   Equtaion  (5)) 

Bl(ITJ)=V(IT,J)»DCOS(THETA(ITJ))»DCOS(PHI(ITJ)) 
B2(IT/])=V(IT,J)*DCOS(THETA(IT,J))*DSIN(PHI(IT,J)) 
B3(IT,J)=V(IT#J)*DSIN(THETA(rrj)) 

70  CONTINUE 

80    CONTINUE 
DO90I=l,LEN 

90    CONTINUE 

Calculate  DELTA  and  TNOT. 


D2=0 

T1=0 

DEN=0 

DO  200  IT=1,LEN 

Dl=(((Bl(rr,l)-Bl(rr,2))*(Bl(IT/l)-Bl(IT,2)))+ 
((B2(IT,3  )-B2(IT,2)),(B2(IT/l)-B2(IT^)))+ 
((B3(rr,l)-B3(rr,2)),(B3(IT,l)-B3(IT,2)))) 

DEN=DEN+D1 

T1=T1+(IPC(IT)»D1) 

D2=D2+(((B1  (IT,1  )-Bl  (IT,2))*(Y1  (IT,!  )-Yl  (IT,2)))+ 

((B2(rr,i)-B2(rr,2))*(Y2(iT/i)-Y2(rr/2)))+ 

(B3(IT,l)-B3(IT,2)m3(IT,l)-Y3(IT,2)))) 

c(rr/i)  =  BiaT,2)-Bi(iT/i) 

C(IT,2)  =  B2(IT,2)  -  62(1X1) 
C(IT,3)  =  B3(IT,2)-B3(IT/1) 

200    CONTINUE 

TD  =  DSQRT(DEN) 
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DEN=DEN/LEN 

TNOT=(Tl/LEN)/DEN 

DELTA=-(D2/LEN)/DEN 


Calculate  m  and  g. 


D1=O.ODO 
D2=O.0D0 
DO210IT=l,LEN 

DD=(<(B1(IT,1)-B1(IT,2))»(B1(IT,1)-BHIT,2)))+ 
(B2(IT,1  )-B2(IT,2))*(B2(IT,l  )-B2(IT,2)))+ 
((B3(IT,1  )-B3(IT,2)HB3(IT,l  )-B3(U ,!))))' 
(IPC(IT)-TNOT)*(IPC(IT)-TNOT)) 
Dl  =  Dl  +  DD 

D2=D2+(((Bl(rr/l)-Bl(IT/2))*(Yl(IT/l)-YHIT,2)))+ 
((B2(IT,1)-B2(IT/2))»(Y2(IT,1)-Y2(IT,2)))+ 
(B3(IT/l)-B3(IT/2))*(Y3(IT/l)-Y3(IT/2))))» 
(IPC(IT)-TNOT) 
210    CONTINUE 

MEQ=-<D2/LEN)/(D1  /LEN) 

TM  =  D1 

DO230J=l,2 

DO220IT=l,LEN 

YC1  0T,J)=Y1  (IT,J)+(B1  ([T  J)*(DELT  A+(MEQ*(IPC(IT)-TNOT)))) 
YC2([T/J)=Y2(rr,J)+(B2(IT,J)*(DELTA+(MEQ»(IPC(rr)-TNOT)))) 
YC3aT/J)=Y3(ITJ)+(B3(IT/])»(DELTA+(MEQ*(IPC(rr)-TNOT)))) 
CEQ(IT)=DELTA+(MEQ»(IPC(IT)-TNOT)) 
220         CONTINUE 
230    CONTINUE 


Create  the  output  files. 


SSR  =  0.0DO 
SSB  =  0.0D0 
DO  231  IT  =  1,LEN 

SSR  =  SSR  +  (YC1(IT,1)  -  YC1(IT,2))"2  +  (YC2(IT,1)  - 

YC2(IT,2))"2  +  (YC3(n\l)  -  YC3(IT,2))**2 
Dl=(((Bl(IT,l)-Bl(IT,2))»(Bl(IT,l)-Bl(IT,2)))+ 
((B2(IT,1  )-B2(IT,2))*(B2(IT,l  )-B2(IT,2))H 
((B3(IT,1)-B3(IT/2))*(B3(IT/1)-B3(IT,2)))) 
SSB  =  SSB  +  DrGEQ(IT)"2 
231     CONTINUE 

MSR  =  SSR/CLEN  -  2) 

MSB  =  SSB/2 

S1G  =  DSQRT(MSR) 

TD  =  TD*DELTA/SIC 

TM  =  MEQ*DSQRT(TM)/S1C 


Calculate  VDEL,SDEL,VM,SDVM    (  variances  and  standard  deviations  ). 
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DUM2  =  O.ODO 
DUM3  =  O.ODO 
DUM5  =  O.ODO 
DUM6  =  O.ODO 
DO  235  1=1,3 
D0  233J  =  1,3 
DUM1  =  O.ODO 
DUM4  =  O.ODO 
D0  232IT  =  1,LEN 

DUM1  =  DUM1  +  (COT,irC(IT,J)) 

DUM4  =  DUM4+((IPC(IT)-TNOT)"2  *C(IT,I)»C(n\J)) 

232  CONTINUE 

DUM2  =  DUM2  +  (M(I,J)  *  DUM1) 
DUM5  =  DUM5  +  (M(IJ)  •  DUM4) 

233  CONTINUE 

DO  234  IT  =  1,LEN 

DUM3  =  DUM3  +  (C(IT,I)*»2) 

DUM6  =  DUM6  +  ((IPC(IT)  -  TNOT)"2  •  C(IT,I)~2) 

234  CONTINUE 

235  CONTINUE 

VDEL  =  DUM2/DUM3"2 
SDEL  =  DSQRT(VDEL) 
VM  =  DUM5/DUM6"2 
SDVM  =  DSQRT(VM) 

DUM2  =  O.ODO 
DO  238  IT  =  1,LEN 
DUM1  =  O.ODO 
DO  237  I  =  1,3 
DO  236  J  =  1,3 

DUM1  =  DUM1  +  (C(IT,I)»C(IT,J)*M(I,J)) 

236  CONTINUE 

237  CONTINUE 

DUM2  =  DUM2  +  ((IPC(IT)  -TNOT)  *  DUM1) 

238  CONTINUE 

COVDM  =  DUM2ADUM3  *  DUM6) 

270    READai,'(A)',END=280)LINE 

GOTO  270 
280    BACKSPACE  11 

IF(SDEL.EQ.O.ODO)  THEN 

DELTAN  =  SDEL 

MEQN  =  SDEL 

R  =  SDEL 

GOTO  285 
ENDIF 

DELTAN  =  DABS(DELTA)/SDEL 
MEQN  =  DABS(MEQ)/SDVM 
R  =  COVDM/(SDEL*SDVM) 
285     DELTA  =  DELTA  *  1000.0DO 
MEQ  =  MEQMOOO.0D0 
SDEL  =  SDEL  *1000.0DO 

67 


SDVM  =  SDVM  *  1000.0DO 
WRITE(n,330)DELTA,SDEL,DELTAN,MEQ,SDVM,MEQN,R,TNOT,LEN 

350     FORMAT(F5.2,2X,F5.4,2X,F6.2^X,F73/2X,F7.6/2X,2(F72,2X)# 
F8.2,2X,I2) 

WRITE(V)'         PROGRAM  COMPLETED  !' 
END 


SUBROUTINE  TCEN(LL,C,VV,A2,P1,P2,ANCLE,TIME, 
V0,V1,DEPTH,IFLAG) 


TGEN  generates  transit  time  and  elevation  angle  at  a  target  if  given  the  horizontal  range,  the  depth  of 
the  sensor  and  the  target,  the  layer  boundaries  and  the  gradients.   Isogradiant  raytracing  is  used. 


Calling  Arguments  are  as  follows: 

LL  -  An  array  containing  the  layer  midpoints. 

G  -  An  array  containing  the  gradients  for  each  layer. 

V  V  -  An  array  containing  the  velocity  at  each  layer. 

A 2  -  The  depth  of  the  sensor  (  positive  down  ). 

PI  -  Range  of  the  target  (  horizontal  down  ). 

P2  -  Depth  of  the  target  (  positive  down  ). 

V0,V1  -  The  values  for  a  straight  line  single  layer  regression  of  depth  vs.  velocity. 

DEPTH  -  An  array  containing  the  depth  of  each  layer. 

Return  arguments  are  as  follows: 

ANGLE      -  The  final  angle  at  the  target. 
TIME  -  The  time  of  transit. 


User  Notes: 

All  floating  point  numbers  are  defines  :  REAL*8 
All  times  are  in  seconds,  and  all  angles  in  radians. 


DIMENSION  L(55),G(55),V(55),LL(55),W(55) 

DIMENSION  TH(55),T(55),VZ(55) 

DIMENSION  C2(55),TT(55),DEPTH(55) 

REAL*8  L,C,LL,VV,TH,VZ,C2,TT,T,R0,A] 

REAL*8  A2,P1,P2,C1,C22,THETA,VM 

REAL*8  THETAZ,RV,R,T1ME,ANGLE,EP,DZ,DEPTH,V0,V1,V 


Initialization:  Set  the  value  for  DZ,  the  layer  thickness.  The  sensor  is  assumed  to  be  at 
RANGE  0.   Determine  the  values  for  J,  which  is  1+Number  of  layers  less  than  or  equal  to  the 
sensor  depth,  and  I,  which  is  the  number  of  layers  less  than  or  equal  to  the  torpedo  depth. 
Redefine  the  endpoints  of  those  layers  locally  to  be  the  depths  of  the  torpedo  and  sensor. 
Define  local  values  for  the  LL  and  W  arravs. 
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EP=0.1D-6 

DZ=25.0D0 

A1=O.ODO 

H 

1=0 
DO10K=l,55 

UK)=LUK) 

V(K)=W(K) 

IF(DEPTH(K).LE.A2)  J=J+1 
IF(DEPTH(K).LE.P2)  1=1+1 
10    CONTINUE 

IF(l.LE.O)  1=1 

N=l+J-1 

V(I)=V(I)+C(I)*(P2-UI)) 

V(J)=V(J-1)+C(J-D*(A2-L(J-1)) 

L(I)=P2 

IF(A2.GT.L0-D)G(J)=G(J-1) 

LG)=A2 

R0=P1-A1 

V1=(V(J)-V(I))/(A2-P2) 

V0=V0)  -  V1*A2 


Calculate  an  initial  estimate  for  the  angle  and  time  using  a  single  layer  approximation. 


C22=-V0/V1 

C1=((0.5D0)»(P1-AD) 

Cl=Cl+((0.5D0)*(L(I)-L(J))*(UI)+La)-2.0D0* 

C22)/(P1-AD) 
IF(IFLAC.CT.O)  GOTO  48 
THETA=DATAN((A1-C1)/(A2-C22)) 
WRITEC,*)'        DEFINE  THETA  AGAIN' 
GOTO  49 

48  THETA=ANGLE 

49  CONTINUE 


Use  the  angle  THETA  to  rayrrace  back  through 
and  the  velocity  to  calculate  the  entrance  andgk 

all  the  layers. 
'  at  each  layer 

First, 

use 

the 

ray 

invariant 

(RV) 

R  =  A1 

Find  the  maximum  sound  speed. 

VM  =  V(I) 
DO  495  K  =  |J 

IF  CV(K).GT.VM)  VM  =  V(K) 
495    CONTINUE 

50     RV=DCOS(THETA)/V(J) 
DO60K=IJ 

TH(K)=DACOS(RV*V(K)) 
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T(K)-DTANfTH(K)) 
VZ(K)=V(K)-L(K)*G(K) 
C2(K)=-VZ(K)/C(K) 
60    CONTINUE 


Using  the  angle  just  calculated,  iterate  backwords  through  the  layers  from  sensor  to  target  to  get 
the  horizontal  range.   Stop  at  the  depth  of  the  target. 

R=0.ODO 

DO  70  K=J,1+1,-1 

Q=R-T(K)»(L(K)-C2(K-1)) 

R=C1  +T(K-1  )*(L(K-1  )-C2(K-l )) 
70    CONTINUE 


Test  if  the  value  for  the  range  is  within  torerance.   If  not,  redefine  THETA,  the  initial  angle, 
and  raytrace  again.   If  vvthin  tolerance,  calculate  the  time  of  travel  based  on  THETA,  and  return. 

EP=0.1D-6 

IF((DABS(R-P1)).LE.EP)  GOTO  100 

THETAZ=THETA 

THETA  =  DATAN(DTAN(THETAZ)*(R-A1)/R0) 

GOTO  50 

100    TT(J)=DLOG((1.0D0+DSIN(TH(I)))/(DCOS(TH(J)))) 
TIME=0.0DO 
DO110K=J-l,I,-l 

TT(K)=DLOG((1.0D0+DSIN(TH(K)))/(DCOS(TH(K)))) 
TIME=TIME+(TT(K)-TT(K+1))/G(K) 
110    CONTINUE 

ANGLE=THETA 

RETURN 

END 
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